orsa_interaction_MPI.cc

Go to the documentation of this file.
00001 /* 
00002    ORSA - Orbit Reconstruction, Simulation and Analysis
00003    Copyright (C) 2002-2004 Pasquale Tricarico
00004    
00005    This program is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU General Public License
00007    as published by the Free Software Foundation; either version 2
00008    of the License, or (at your option) any later version.
00009    
00010    As a special exception, Pasquale Tricarico gives permission to
00011    link this program with Qt commercial edition, and distribute the
00012    resulting executable, without including the source code for the Qt
00013    commercial edition in the source distribution.
00014    
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019    
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the Free Software
00022    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00023 */
00024 
00025 #include "orsa_interaction.h"
00026 #include "orsa_secure_math.h"
00027 #include "orsa_universe.h"
00028 
00029 #include <cmath>
00030 #include <iostream>
00031 
00032 #include <mpi.h>
00033 
00034 using namespace std;
00035 
00036 namespace orsa {
00037   
00038   // Newton_MPI
00039   
00040   Newton_MPI::Newton_MPI() : Interaction() {
00041     av = av_local = bv = 0;
00042   }
00043   
00044   Newton_MPI::Newton_MPI(const Newton_MPI &) : Interaction() {
00045     av = av_local = bv = 0;
00046   }
00047   
00048   Newton_MPI::~Newton_MPI() {
00049     free(bv);
00050     free(av);
00051     free(av_local);
00052   }
00053   
00054   Interaction * Newton_MPI::clone() const {
00055     return new Newton_MPI(*this);
00056   }
00057   
00058   void Newton_MPI::Acceleration(const Frame &fr, vector<Vector> &a) {
00059     
00060     int rank, size;
00061     MPI_Comm_rank(MPI_COMM_WORLD, &rank );
00062     MPI_Comm_size(MPI_COMM_WORLD, &size );
00063     
00064     // fprintf(stderr,"inside Newton_MPI::Acceleration() ---> rank = %i\n",rank);
00065     
00066     unsigned int i,j;
00067     
00068     Vector d;
00069     
00070     double ls;
00071     
00072     // sync
00073     MPI_Barrier(MPI_COMM_WORLD); 
00074     
00075     unsigned int frame_size=0;
00076     if (rank==0) {
00077       frame_size = fr.size();
00078     }
00079     //
00080     MPI_Bcast(&frame_size, 1, MPI_INT, 0, MPI_COMM_WORLD);
00081     
00082     if (rank==0) {
00083       g = GetG();
00084     }
00085     //
00086     MPI_Bcast(&g, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00087     
00088     // cerr << "rank: " << rank << "  frame size: " << frame_size << endl;
00089     
00090     if (frame_size < 2) return;
00091     
00092     a.resize(frame_size);
00093     for (i=0;i<frame_size;++i)
00094       a[i].Set(0.0,0.0,0.0);
00095     
00096     // debug
00097     /* {
00098        cerr << "a reset:" << endl;
00099        unsigned int k=0;
00100        while (k<frame_size) {
00101        fprintf(stderr,"a[%i].Length() = %g\n",k,a[k].Length());
00102        ++k;
00103        }
00104        }
00105     */
00106     
00107     // 4 = mass (1) + position (3)
00108     bv = (double *)realloc(bv,4*frame_size*sizeof(double)); 
00109     
00110     if (rank==0) {
00111       unsigned int k=0;
00112       while (k<frame_size) {
00113         bv[4*k]   = fr[k].mass();
00114         //
00115         bv[4*k+1] = fr[k].position().x;
00116         bv[4*k+2] = fr[k].position().y;
00117         bv[4*k+3] = fr[k].position().z;
00118         //
00119         ++k;
00120       }
00121     }
00122     
00123     // debug
00124     /* {
00125        cerr << "initialized bv:" << endl;
00126        unsigned int k=0;
00127        while (k<4*frame_size) {
00128        fprintf(stderr,"bv[%i] = %g\n",k,bv[k]);
00129        ++k;
00130        }
00131        }
00132     */
00133     
00134     // cerr << "before Bcast..." << endl;
00135     
00136     MPI_Bcast(bv, 4*frame_size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00137     
00138     // cerr << "after Bcast..." << endl;
00139     
00140     // debug
00141     /* {
00142        cerr << "bv after Bcast:" << endl;
00143        unsigned int k=0;
00144        while (k<4*frame_size) {
00145        fprintf(stderr,"bv[%i] = %g\n",k,bv[k]);
00146        ++k;
00147        }
00148        }
00149     */
00150     
00151     // print(ff);
00152     
00153     for (i=1+rank;i<frame_size;i+=size) {
00154       
00155       for (j=0;j<i;++j) {
00156         
00157         // fprintf(stderr,"rank: %i   computing %i <--> %i\n",rank,i,j); 
00158         
00159         // if ((ff[i].mass()==0) && (ff[j].mass()==0)) continue;
00160         if ((bv[4*i]==0.0) && (bv[4*j]==0.0)) continue;
00161         
00162         // d = ff[i].DistanceVector(ff[j]);
00163         d.Set(bv[4*j+1]-bv[4*i+1],
00164               bv[4*j+2]-bv[4*i+2],
00165               bv[4*j+3]-bv[4*i+3]);
00166         
00167         // fprintf(stderr,"rank: %i   computing %i <--> %i  d.Length(): %g\n",rank,i,j,d.Length()); 
00168         
00169         ls = d.LengthSquared();
00170         
00171         // fprintf(stderr,"rank: %i   computing %i <--> %i  ls: %g\n",rank,i,j,ls); 
00172         
00173         if (d.IsZero()) {
00174           cerr << "*** Warning: two objects in the same position!" << endl;
00175           continue;
00176         }
00177         
00178         d *= secure_pow(ls,-1.5);
00179         
00180         // fprintf(stderr,"rank: %i   computing %i <--> %i  modified-d.Length(): %g\n",rank,i,j,d.Length()); 
00181         
00182         // fprintf(stderr,"rank: %i   computing %i <--> %i  pre-change a[%i].Length(): %g   a[%i].Length(): %g\n",rank,i,j,i,a[i].Length(),j,a[j].Length()); 
00183         
00184         // a[i] += d * ff[j].mass();
00185         // a[j] -= d * ff[i].mass();
00186         // a[i] += d * bv[4*j];
00187         // a[j] -= d * bv[4*i];
00188         a[i] = a[i] + d * bv[4*j];
00189         a[j] = a[j] - d * bv[4*i];
00190         
00191         // fprintf(stderr,"rank: %i   computing %i <--> %i  end of cycle a[%i].Length(): %g   a[%i].Length(): %g\n",rank,i,j,i,a[i].Length(),j,a[j].Length()); 
00192       }  
00193     } 
00194     
00195     // cerr << "rank " << rank << "    out of the loop..." << endl;
00196     
00197     // double *av = (double *)malloc(3*frame_size*sizeof(double));
00198     av = (double *)realloc(av,3*frame_size*sizeof(double));
00199     
00200     if (rank == 0) {
00201       unsigned int k=0;
00202       while (k<3*frame_size) {
00203         av[k] = 0.0;
00204         ++k;
00205       }
00206     }
00207     
00208     // double *av_local = (double *)malloc(3*frame_size*sizeof(double));
00209     av_local = (double *)realloc(av_local,3*frame_size*sizeof(double));
00210     
00211     {
00212       unsigned int k=0;
00213       while (k<frame_size) {
00214         av_local[3*k]   = a[k].x;
00215         av_local[3*k+1] = a[k].y;
00216         av_local[3*k+2] = a[k].z;
00217         ++k;
00218       }
00219     }
00220     
00221     // debug   
00222     /* {
00223        unsigned int k=0;
00224        while (k<3*frame_size) {
00225        fprintf(stderr,"rank: %i  av_local[%i]: %g\n",rank,k,av_local[k]);
00226        ++k;
00227        }
00228        }
00229     */
00230     
00231     MPI_Reduce(av_local, av, 3*frame_size, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00232     
00233     if (rank == 0) {
00234       unsigned int k=0;
00235       while (k<frame_size) {
00236         a[k].x = av[3*k];
00237         a[k].y = av[3*k+1];
00238         a[k].z = av[3*k+2];     
00239         //
00240         a[k] *= g;
00241         //
00242         // fprintf(stderr,"acc[%i].Length() = %g\n",k,a[k].Length());
00243         //
00244         ++k;
00245       }
00246     }
00247   }
00248   
00249   double Newton_MPI::PotentialEnergy(const Frame &f) {
00250     Newton newton;
00251     return (newton.PotentialEnergy(f));
00252   }
00253   
00254 } // namespace orsa

Generated on Thu Jul 13 06:45:22 2006 for liborsa by  doxygen 1.4.7