ExaFMM 1
Fast-multipole Method for exascale systems
|
00001 /* 00002 Copyright (C) 2011 by Rio Yokota, Simon Layton, Lorena Barba 00003 00004 Permission is hereby granted, free of charge, to any person obtaining a copy 00005 of this software and associated documentation files (the "Software"), to deal 00006 in the Software without restriction, including without limitation the rights 00007 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 00008 copies of the Software, and to permit persons to whom the Software is 00009 furnished to do so, subject to the following conditions: 00010 00011 The above copyright notice and this permission notice shall be included in 00012 all copies or substantial portions of the Software. 00013 00014 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 00015 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 00016 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 00017 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 00018 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 00019 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 00020 THE SOFTWARE. 00021 */ 00022 #if VTK 00023 #define VTK_EXCLUDE_STRSTREAM_HEADERS 00024 #include <vtkPoints.h> 00025 #include <vtkGraphLayoutView.h> 00026 #include <vtkMutableUndirectedGraph.h> 00027 #include <vtkRenderWindow.h> 00028 #include <vtkRenderWindowInteractor.h> 00029 #include <vtkInteractorStyleTrackballCamera.h> 00030 #include <vtkGraphWriter.h> 00031 #include <vtkViewTheme.h> 00032 #endif 00033 #include "parallelfmm.h" 00034 00035 struct JVertex { 00036 int Ista; 00037 int Iend; 00038 vect X; 00039 }; 00040 00041 struct Vertex : JVertex { 00042 #if VTK 00043 vtkIdType Id; 00044 #endif 00045 vect F; 00046 }; 00047 00048 const std::string INPUT_PATH = "../../fdgl_in/"; 00049 const std::string OUTPUT_PATH = "../../fdgl_out/"; 00050 00051 std::vector<int> edges; 00052 std::vector<Vertex> vertices; 00053 std::vector<Vertex> verticesG; 00054 std::vector<JVertex> jvertices; 00055 typedef std::vector<Vertex>::iterator V_iter; 00056 typedef std::vector<JVertex>::iterator JV_iter; 00057 #if VTK 00058 vtkMutableUndirectedGraph *graph = vtkMutableUndirectedGraph::New(); 00059 vtkGraphLayoutView *view = vtkGraphLayoutView::New(); 00060 vtkGraphWriter *writer = vtkGraphWriter::New(); 00061 #endif 00062 int numVertices, numEdges, localVertices, maxEdges=0; 00063 int *offset = new int [MPISIZE]; 00064 static real radius = 2; 00065 00066 double get_time() { 00067 struct timeval tv; 00068 gettimeofday(&tv,NULL); 00069 return double(tv.tv_sec+tv.tv_usec*1e-6); 00070 } 00071 00072 void splitRange(int &begin, int &end, int iSplit, int numSplit) { 00073 int size = end - begin; 00074 int increment = size / numSplit; 00075 int remainder = size % numSplit; 00076 begin += iSplit * increment + std::min(iSplit,remainder); 00077 end = begin + increment; 00078 if( remainder > iSplit ) end++; 00079 } 00080 00081 void initVertices() { 00082 int begin = 0; 00083 int end = numVertices; 00084 splitRange(begin,end,MPIRANK,MPISIZE); 00085 localVertices = end - begin; 00086 srand48(0); 00087 for( int i=0; i<numVertices; ++i ) { 00088 Vertex vertex; 00089 #if VTK 00090 vertex.Id = graph->AddVertex(); 00091 #endif 00092 vertex.Ista = vertex.Iend = 0; 00093 vertex.X[0] = 2 * drand48() - 1; 00094 vertex.X[1] = 2 * drand48() - 1; 00095 vertex.X[2] = 2 * drand48() - 1; 00096 vertex.F = 0; 00097 if( begin <= i && i < end ) vertices.push_back(vertex); 00098 } 00099 verticesG.resize(numVertices); 00100 } 00101 00102 void readEdges(std::ifstream &fid) { 00103 std::string line0, line1; 00104 int begin = 0; 00105 int end = numVertices; 00106 splitRange(begin,end,MPIRANK,MPISIZE); 00107 fid >> line0 >> line1; 00108 while( atoi(line0.c_str()) < begin ) { 00109 fid >> line0 >> line1; 00110 } 00111 int prevSrc=0, crntSrc=-1; 00112 while( !fid.eof() && atoi(line0.c_str()) < end ) { 00113 crntSrc = atoi(line0.c_str()) - begin; 00114 if (crntSrc != prevSrc) { 00115 vertices[prevSrc].Iend = edges.size(); 00116 maxEdges = std::max(vertices[prevSrc].Iend - vertices[prevSrc].Ista, maxEdges); 00117 vertices[crntSrc].Ista = edges.size(); 00118 prevSrc = crntSrc; 00119 } 00120 edges.push_back(atoi(line1.c_str())); 00121 fid >> line0 >> line1; 00122 } 00123 vertices[prevSrc].Iend = edges.size(); 00124 maxEdges = std::max(vertices[prevSrc].Iend - vertices[prevSrc].Ista, maxEdges); 00125 } 00126 00127 void readGraph(std::string fileid) { 00128 std::string line0, line1; 00129 std::string file=INPUT_PATH+fileid+".graph"; 00130 std::ifstream fid(file.c_str(),std::ios::in); 00131 // First line in file contains <num nodes> <num edges> 00132 fid >> line0 >> line1; 00133 numVertices = atoi(line0.c_str()); 00134 numEdges = atoi(line1.c_str()); 00135 if (numVertices > 0) { // handling bad input files 00136 initVertices(); 00137 readEdges(fid); 00138 } else { 00139 std::cerr << "Bad input file: " << file << std::endl; 00140 } 00141 jvertices.resize(numEdges); 00142 fid.close(); 00143 } 00144 00145 void shiftVertices() { 00146 int newSize = 0; 00147 int oldSize = vertices.size(); 00148 const int bytes = sizeof(Vertex) / 8; 00149 const int isend = (MPIRANK + 1 ) % MPISIZE; 00150 const int irecv = (MPIRANK - 1 + MPISIZE) % MPISIZE; 00151 MPI_Request sreq,rreq; 00152 00153 MPI_Isend(&oldSize,1,MPI_INT,irecv,0,MPI_COMM_WORLD,&sreq); 00154 MPI_Irecv(&newSize,1,MPI_INT,isend,0,MPI_COMM_WORLD,&rreq); 00155 MPI_Wait(&sreq,MPI_STATUS_IGNORE); 00156 MPI_Wait(&rreq,MPI_STATUS_IGNORE); 00157 00158 std::vector<Vertex> buffer(newSize); 00159 MPI_Isend(&vertices[0],oldSize*bytes,MPI_DOUBLE,irecv, 00160 1,MPI_COMM_WORLD,&sreq); 00161 MPI_Irecv(&buffer[0] ,newSize*bytes,MPI_DOUBLE,isend, 00162 1,MPI_COMM_WORLD,&rreq); 00163 MPI_Wait(&sreq,MPI_STATUS_IGNORE); 00164 MPI_Wait(&rreq,MPI_STATUS_IGNORE); 00165 vertices = buffer; 00166 } 00167 00168 void shiftEdges() { 00169 int newSize = 0; 00170 int oldSize = edges.size(); 00171 const int isend = (MPIRANK + 1 ) % MPISIZE; 00172 const int irecv = (MPIRANK - 1 + MPISIZE) % MPISIZE; 00173 MPI_Request sreq,rreq; 00174 00175 MPI_Isend(&oldSize,1,MPI_INT,irecv,0,MPI_COMM_WORLD,&sreq); 00176 MPI_Irecv(&newSize,1,MPI_INT,isend,0,MPI_COMM_WORLD,&rreq); 00177 MPI_Wait(&sreq,MPI_STATUS_IGNORE); 00178 MPI_Wait(&rreq,MPI_STATUS_IGNORE); 00179 00180 std::vector<int> buffer(newSize); 00181 MPI_Isend(&edges[0] ,oldSize,MPI_INT,irecv, 00182 1,MPI_COMM_WORLD,&sreq); 00183 MPI_Irecv(&buffer[0],newSize,MPI_INT,isend, 00184 1,MPI_COMM_WORLD,&rreq); 00185 MPI_Wait(&sreq,MPI_STATUS_IGNORE); 00186 MPI_Wait(&rreq,MPI_STATUS_IGNORE); 00187 edges = buffer; 00188 } 00189 00190 void gatherVertices() { 00191 int sendCnt = localVertices * sizeof(Vertex) / 8; 00192 int *recvCnt = new int [MPISIZE]; 00193 MPI_Allgather(&sendCnt,1,MPI_INT,recvCnt,1,MPI_INT,MPI_COMM_WORLD); 00194 offset[0] = 0; 00195 for( int i=1; i!=MPISIZE; ++i ) { 00196 offset[i] = offset[i-1] + recvCnt[i-1]; 00197 } 00198 MPI_Allgatherv(&vertices[0],sendCnt,MPI_DOUBLE,&verticesG[0],recvCnt,offset,MPI_DOUBLE,MPI_COMM_WORLD); 00199 delete[] recvCnt; 00200 } 00201 00202 #if VTK 00203 void setEdges() { 00204 gatherVertices(); 00205 for( int irank=0; irank!=MPISIZE; ++irank ) { 00206 for( V_iter V=vertices.begin(); V!=vertices.end(); ++V ) { 00207 for( int i=V->Ista; i<V->Iend; ++i ) { 00208 // Data comes in duplicates - only add src <= dest 00209 if (V->Id < verticesG[edges[i]].Id) { 00210 graph->AddEdge(V->Id, verticesG[edges[i]].Id); 00211 } 00212 } 00213 } 00214 shiftVertices(); 00215 shiftEdges(); 00216 } 00217 } 00218 00219 void setVertices() { 00220 vtkPoints *points = vtkPoints::New(); 00221 for( int irank=0; irank!=MPISIZE; ++irank ) { 00222 for( V_iter V=vertices.begin(); V!=vertices.end(); ++V ) { 00223 points->InsertNextPoint(V->X[0], V->X[1], V->X[2]); 00224 } 00225 shiftVertices(); 00226 } 00227 graph->SetPoints(points); 00228 } 00229 #endif 00230 00231 void repulsion(ParallelFMM<Laplace> &FMM) { 00232 Bodies bodies(localVertices); 00233 Bodies jbodies; 00234 Cells cells, jcells; 00235 B_iter B = bodies.begin(); 00236 vect Xmin = B->X; 00237 vect Xmax = B->X; 00238 for( V_iter V=vertices.begin(); V!=vertices.end(); ++V, ++B ) { 00239 B->X = V->X; 00240 B->SRC = -100; 00241 B->TRG = 0; 00242 B->IBODY = V-vertices.begin(); 00243 B->IPROC = MPIRANK; 00244 B->ICELL = 0; 00245 if( Xmin[0] < B->X[0] ) Xmin[0] = B->X[0]; 00246 if( Xmin[1] < B->X[1] ) Xmin[1] = B->X[1]; 00247 if( Xmin[2] < B->X[2] ) Xmin[2] = B->X[2]; 00248 if( Xmax[0] > B->X[0] ) Xmax[0] = B->X[0]; 00249 if( Xmax[1] > B->X[1] ) Xmax[1] = B->X[1]; 00250 if( Xmax[2] > B->X[2] ) Xmax[2] = B->X[2]; 00251 } 00252 radius = std::sqrt(norm(Xmax-Xmin)); 00253 FMM.setGlobDomain(bodies); 00254 FMM.octsection(bodies); 00255 FMM.bottomup(bodies,cells); 00256 FMM.commBodies(cells); 00257 jbodies = bodies; 00258 jcells = cells; 00259 FMM.commCells(jbodies,jcells); 00260 FMM.downward(cells,jcells); 00261 FMM.unpartition(bodies); 00262 for( B=bodies.begin(); B!=bodies.end(); ++B ) { 00263 vertices[B->IBODY].F[0] = B->TRG[1]; 00264 vertices[B->IBODY].F[1] = B->TRG[2]; 00265 vertices[B->IBODY].F[2] = B->TRG[3]; 00266 } 00267 } 00268 00269 void commVertices() { 00270 gatherVertices(); 00271 for( V_iter VI=vertices.begin(); VI!=vertices.end(); ++VI ) { 00272 for( int i=VI->Ista; i<VI->Iend; ++i ) { 00273 V_iter VJ = verticesG.begin()+edges[i]; 00274 jvertices[i].Ista = VJ->Ista; 00275 jvertices[i].Iend = VJ->Iend; 00276 jvertices[i].X = VJ->X; 00277 } 00278 } 00279 } 00280 00281 void spring() { 00282 float l = 2 / pow(numVertices,1./3); 00283 for( V_iter VI=vertices.begin(); VI!=vertices.end(); ++VI ) { 00284 for( int i=VI->Ista; i<VI->Iend; ++i ) { 00285 JV_iter VJ = jvertices.begin()+i; 00286 vect dist = VI->X - VJ->X; 00287 float R = sqrtf(norm(dist) + EPS2); 00288 float weight = (VI->Iend-VI->Ista) * (VJ->Iend-VJ->Ista); 00289 VI->F -= dist / R * (R - l / weight); 00290 } 00291 } 00292 } 00293 00294 void moveVertices() { 00295 for( V_iter V=vertices.begin(); V!=vertices.end(); ++V ) { 00296 vect dX; 00297 if( norm(V->F) < EPS ) dX = 0; 00298 else dX = V->F / std::sqrt(norm(V->F)); 00299 V->X += dX; 00300 } 00301 } 00302 00303 #if VTK 00304 void drawGraph() { 00305 vtkViewTheme* theme = vtkViewTheme::New(); 00306 theme->SetBackgroundColor(1.,1.,1.); 00307 theme->SetBackgroundColor2(1.,1.,1.); 00308 theme->SetPointColor(.2,.2,.2); 00309 theme->SetCellColor(.2,.2,.2); 00310 view->ApplyViewTheme(theme); 00311 theme->Delete(); 00312 view->AddRepresentationFromInput(graph); 00313 view->SetLayoutStrategy("Pass Through"); 00314 view->GetInteractor()->GetRenderWindow()->SetSize(700,700); 00315 view->ResetCamera(); 00316 view->Render(); 00317 } 00318 00319 void finalizeGraph() { 00320 vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New(); 00321 view->GetInteractor()->SetInteractorStyle(style); 00322 view->GetInteractor()->Initialize(); 00323 view->GetInteractor()->Start(); 00324 } 00325 00326 void writeGraph(std::string fileid, int step) { 00327 std::stringstream ss; 00328 ss << OUTPUT_PATH << fileid << '_' << step << ".vtk"; 00329 writer->SetFileName(ss.str().c_str()); 00330 writer->SetInput(graph); 00331 writer->Write(); 00332 } 00333 #endif 00334 00335 int main() { 00336 double t0, t[7] = {0,0,0,0,0,0,0}; 00337 IMAGES = 0; 00338 THETA = 0.6; 00339 ParallelFMM<Laplace> FMM; 00340 FMM.initialize(); 00341 std::string fileid="jagmesh1_v936_e2664"; 00342 // std::string fileid="add32_v4960_e9462"; 00343 // std::string fileid="4elt_v15606_e45878"; 00344 // std::string fileid="finan512_v74752_e261120"; 00345 // std::string fileid="uni-socialgraph-anonymized_v59216214_e92522012"; 00346 // std::string fileid="mhrw-socialgraph-anonymized_v72261577_e159189465"; 00347 // std::string fileid="uplrg_v1000_e2000"; 00348 // std::string fileid="uplrg_v10000_e20000"; 00349 // std::string fileid="uplrg_v100000_e200000"; 00350 // std::string fileid="uplrg_v1000000_e2000000"; 00351 // std::string fileid="uplrg_v10000000_e20000000"; 00352 // std::string fileid="uplrg_v100000000_e200000000"; 00353 00354 t0 = get_time(); 00355 readGraph(fileid); 00356 #if VTK 00357 setEdges(); 00358 #endif 00359 t[0] += get_time() - t0; 00360 00361 for( int step=0; step<1000; ++step ) { 00362 if( MPIRANK == 0 ) std::cout << step << std::endl; 00363 00364 t0 = get_time(); 00365 repulsion(FMM); 00366 t[1] += get_time() - t0; 00367 00368 t0 = get_time(); 00369 commVertices(); 00370 spring(); 00371 t[2] += get_time() - t0; 00372 00373 t0 = get_time(); 00374 moveVertices(); 00375 t[3] += get_time() - t0; 00376 00377 #if VTK 00378 t0 = get_time(); 00379 setVertices(); 00380 t[4] += get_time() - t0; 00381 00382 t0 = get_time(); 00383 if( MPIRANK == 0 && step % 10 == 0 ) drawGraph(); 00384 t[5] += get_time() - t0; 00385 00386 t0 = get_time(); 00387 // if( MPIRANK == 0 ) writeGraph(fileid,step); 00388 t[6] += get_time() - t0; 00389 #endif 00390 00391 } 00392 00393 if( MPIRANK == 0 ) { 00394 std::cout << "V : " << numVertices << std::endl; 00395 std::cout << "E : " << numEdges / 2 << std::endl; 00396 std::cout << "initialize : " << t[0] << std::endl; 00397 std::cout << "repulsion : " << t[1] << std::endl; 00398 std::cout << "spring : " << t[2] << std::endl; 00399 std::cout << "move : " << t[3] << std::endl; 00400 #if VTK 00401 std::cout << "setVer : " << t[4] << std::endl; 00402 std::cout << "draw : " << t[5] << std::endl; 00403 std::cout << "writeGraph : " << t[6] << std::endl; 00404 #endif 00405 } 00406 00407 #if VTK 00408 if( MPIRANK == 0 ) finalizeGraph(); 00409 #endif 00410 FMM.finalize(); 00411 }