1 #ifndef ALUGRID_ZOLTAN_H_INCLUDED 2 #define ALUGRID_ZOLTAN_H_INCLUDED 3 4 #include <iostream> 5 #include <cmath> 6 #include <dune/alugrid/common/alugrid_assert.hh> 7 #include <sstream> 8 9 #include "mpAccess_MPI.h" 10 11 // Warning: Zoltan defines HAVE_MPI and HAVE_PARMETIS itself. However, their definition will 12 // not match ours. The following complicated preprocessor code tries 13 // to cope with this problem. 14 15 #if HAVE_ZOLTAN 16 17 // if DUNE was built with MPI 18 #if HAVE_MPI 19 // undefine our definition of HAVE_MPI before including zoltan_cpp.h 20 #undef HAVE_MPI 21 #define HAVE_MPI_WAS_UNDEFED_HERE 22 #endif 23 24 #if HAVE_PARMETIS 25 #undef HAVE_PARMETIS 26 #define HAVE_PARMETIS_WAS_UNDEFED_HERE 27 #endif 28 29 // include Zoltan's C++ header 30 #include <zoltan_cpp.h> 31 32 // undefine any definition of HAVE_MPI made by Zoltan 33 #ifdef HAVE_MPI 34 #undef HAVE_MPI 35 #endif // #ifdef HAVE_MPI 36 37 // undefine any definition of HAVE_PARMETIS made by Zoltan 38 #ifdef HAVE_PARMETIS 39 #undef HAVE_PARMETIS 40 #endif 41 42 #ifdef HAVE_MPI_WAS_UNDEFED_HERE 43 #ifdef ENABLE_MPI 44 // redefine our definition of HAVE_MPI if it was undef'd before 45 #define HAVE_MPI ENABLE_MPI 46 #else 47 #define HAVE_MPI 1 48 #endif 49 #undef HAVE_MPI_WAS_UNDEFED_HERE 50 #endif // #ifdef HAVE_MPI_WAS_UNDEFED_HERE 51 52 #ifdef HAVE_PARMETIS_WAS_UNDEFED_HERE 53 // redefine our definition of HAVE_PARMETIS if it was undef'd before 54 #define HAVE_PARMETIS ENABLE_PARMETIS 55 #undef HAVE_PARMETIS_WAS_UNDEFED_HERE 56 #endif // #ifdef HAVE_PARMETIS_WAS_UNDEFED_HERE 57 58 #endif // #if HAVE_ZOLTAN 59 60 namespace ALUGridZoltan 61 { 62 using ::ALUGrid::MpAccessGlobal; 63 using ::ALUGrid::MpAccessLocal; 64 65 #if HAVE_ZOLTAN 66 template < class ldb_vertex_map_t, class ldb_edge_set_t > 67 class ObjectCollection 68 { 69 int _rank; 70 ldb_vertex_map_t& _vertexMap ; 71 ldb_edge_set_t& _edgeMap ; 72 typedef typename ldb_edge_set_t::const_iterator edgeType; 73 std::vector< std::vector< std::pair<edgeType,bool> > > _edges; 74 static const int dimension = 3 ; 75 76 public: 77 // constructor ObjectCollection(int rank,ldb_vertex_map_t & vertexMap,ldb_edge_set_t & edgeMap)78 ObjectCollection( int rank, ldb_vertex_map_t& vertexMap, ldb_edge_set_t& edgeMap ) 79 : _rank( rank ), 80 _vertexMap( vertexMap ), 81 _edgeMap( edgeMap ), 82 _edges(0) 83 {} 84 rank()85 int rank() { return _rank; } vertexMap()86 ldb_vertex_map_t& vertexMap() { return _vertexMap; } edgeMap()87 ldb_edge_set_t& edgeMap() { return _edgeMap; } edges()88 std::vector< std::vector<std::pair<edgeType,bool> > >& edges() { return _edges; } edgeIdx(int i,int k)89 int edgeIdx(int i,int k) { alugrid_assert ( i < (int)_edges.size() && k < (int)_edges[i].size() ); 90 return ( (_edges[i][k].second) ? 91 _edges[i][k].first->rightNode() : 92 _edges[i][k].first->leftNode() ) ; } edgeMaster(int i,int k)93 int edgeMaster(int i,int k) { alugrid_assert ( i < (int)_edges.size() && k < (int)_edges[i].size() ); 94 return ( (_edges[i][k].second) ? 95 _edges[i][k].first->rightMaster() : 96 _edges[i][k].first->leftMaster() ) ; } edgeWeight(int i,int k)97 int edgeWeight(int i,int k) { alugrid_assert ( i < (int)_edges.size() && k < (int)_edges[i].size() ); 98 return _edges[i][k].first->weight(); } 99 100 // query functions that respond to requests from Zoltan 101 get_number_of_objects(void * data,int * ierr)102 static int get_number_of_objects(void *data, int *ierr) 103 { 104 ObjectCollection *objs = static_cast<ObjectCollection *> (data); 105 alugrid_assert ( objs ); 106 *ierr = ZOLTAN_OK; 107 return objs->vertexMap().size(); 108 } 109 get_object_list(void * data,int sizeGID,int sizeLID,ZOLTAN_ID_PTR globalID,ZOLTAN_ID_PTR localID,int wgt_dim,float * obj_wgts,int * ierr)110 static void get_object_list(void *data, int sizeGID, int sizeLID, 111 ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, 112 int wgt_dim, float *obj_wgts, int *ierr) 113 { 114 ObjectCollection *objs = static_cast<ObjectCollection *> (data); 115 alugrid_assert ( objs ); 116 alugrid_assert (wgt_dim==1); 117 *ierr = ZOLTAN_OK; 118 119 ldb_vertex_map_t& vertexMap = objs->vertexMap(); 120 121 // In this example, return the IDs of our objects, 122 int i = 0; 123 typedef typename ldb_vertex_map_t :: iterator iterator ; 124 const iterator end = vertexMap.end(); 125 for ( iterator it = vertexMap.begin(); it != end; ++ it, ++i ) 126 { 127 // std::cout << "[" << objs->rank() << "]: "; 128 // std::cout << "vertex(" << i << ") = " << (*it).first.index() << std::endl; 129 globalID[ i ] = (*it).first.index() ; 130 localID [ i ] = i; 131 if (wgt_dim == 1) 132 obj_wgts[ i ] = (*it).first.weight(); 133 } 134 } 135 get_num_edges_list(void * data,int sizeGID,int sizeLID,int num_obj,ZOLTAN_ID_PTR globalID,ZOLTAN_ID_PTR localID,int * numEdges,int * ierr)136 static void get_num_edges_list(void *data, int sizeGID, int sizeLID, 137 int num_obj, 138 ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, 139 int *numEdges, int *ierr) 140 { 141 ObjectCollection *objs = static_cast<ObjectCollection *> (data); 142 alugrid_assert ( num_obj == (int)objs->vertexMap().size() ); 143 if (num_obj == 0) 144 { 145 *ierr = ZOLTAN_OK; 146 return; 147 } 148 149 for (int i=0;i<num_obj;++i) 150 { 151 numEdges[i] = 0; 152 } 153 objs->edges().resize(num_obj); 154 155 const int rank = objs->rank(); 156 157 typename ldb_edge_set_t::const_iterator iEnd = objs->edgeMap().end(); 158 for (typename ldb_edge_set_t::const_iterator it = objs->edgeMap().begin (); it != iEnd; ++it ) 159 { 160 if ( it->leftMaster() == rank ) 161 { 162 int node = it->leftNode(); 163 int i=0; 164 // std::cout << "[" << objs->rank() << "]: "; 165 // std::cout << "edge(" << node << ") = " << it->rightNode() << " " << it->rightMaster() << std::endl; 166 for (;i<num_obj;++i) 167 if ((int)globalID[i] == node) break; 168 alugrid_assert ( i<num_obj ); 169 alugrid_assert (it->rightMaster() >= 0); 170 alugrid_assert (it->weight() >= 0); 171 ++numEdges[i]; 172 objs->edges()[i].push_back( std::make_pair(it,true) ); 173 } 174 if ( it->rightMaster() == rank ) 175 { 176 int node = it->rightNode(); 177 // std::cout << "[" << objs->rank() << "]: "; 178 // std::cout << "edge(" << node << ") = " << it->leftNode() << " " << it->leftMaster() << std::endl; 179 int i=0; 180 for (;i<num_obj;++i) 181 if ((int)globalID[i] == node) break; 182 alugrid_assert ( i<num_obj ); 183 alugrid_assert (it->leftMaster() >= 0); 184 alugrid_assert (it->weight() >= 0); 185 ++numEdges[i]; 186 objs->edges()[i].push_back( std::make_pair(it,false) ); 187 } 188 } 189 *ierr = ZOLTAN_OK; 190 } get_edge_list(void * data,int sizeGID,int sizeLID,int num_obj,ZOLTAN_ID_PTR globalID,ZOLTAN_ID_PTR localID,int * num_edges,ZOLTAN_ID_PTR nborGID,int * nborProc,int wgt_dim,float * ewgts,int * ierr)191 static void get_edge_list(void *data, int sizeGID, int sizeLID, 192 int num_obj, ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, 193 int *num_edges, 194 ZOLTAN_ID_PTR nborGID, int *nborProc, 195 int wgt_dim, float *ewgts, int *ierr) 196 { 197 ObjectCollection *objs = static_cast<ObjectCollection *> (data); 198 if (num_obj == 0) 199 { 200 *ierr = ZOLTAN_OK; 201 return; 202 } 203 int k=0; 204 for (int j=0;j<num_obj;++j) 205 { 206 for (int l=0;l<num_edges[j];++l) 207 { 208 // std::cout << "[" << objs->rank() << "]: "; 209 // std::cout << "v(" << j << ")(" << l << ")=" << objs->edgeIdx(j,l) << " " << objs->edgeMaster(j,l) << std::endl; 210 nborGID[k] = objs->edgeIdx(j,l); 211 nborProc[k] = objs->edgeMaster(j,l); 212 alugrid_assert ( nborProc[k] >= 0 ); 213 if (wgt_dim==1) 214 { 215 ewgts[k]=objs->edgeWeight(j,l); 216 alugrid_assert ( ewgts[k] >= 0 ); 217 } 218 ++k; 219 } 220 } 221 *ierr = ZOLTAN_OK; 222 } 223 224 // return dimension of coordinates get_num_geometry(void * data,int * ierr)225 static int get_num_geometry(void *data, int *ierr) 226 { 227 *ierr = ZOLTAN_OK; 228 return dimension; 229 } 230 get_geometry_list(void * data,int sizeGID,int sizeLID,int num_obj,ZOLTAN_ID_PTR globalID,ZOLTAN_ID_PTR localID,int num_dim,double * geom_vec,int * ierr)231 static void get_geometry_list(void *data, int sizeGID, int sizeLID, 232 int num_obj, 233 ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, 234 int num_dim, double *geom_vec, int *ierr) 235 { 236 ObjectCollection *objs = static_cast<ObjectCollection *> (data); 237 alugrid_assert ( objs ); 238 239 if ( (sizeGID != 1) || (sizeLID != 1) || (num_dim != dimension)) 240 { 241 *ierr = ZOLTAN_FATAL; 242 return; 243 } 244 245 *ierr = ZOLTAN_OK; 246 247 ldb_vertex_map_t& vertexMap = objs->vertexMap(); 248 249 int idx = 0; 250 typedef typename ldb_vertex_map_t :: iterator iterator ; 251 double coord[ 3 ]; 252 const iterator end = vertexMap.end(); 253 for ( iterator it = vertexMap.begin(); it != end; ++ it ) 254 { 255 (*it).first.computeBaryCenter( coord ); 256 for( int d=0; d<dimension; ++d, ++idx ) 257 geom_vec[ idx ] = coord[ d ]; 258 } 259 } 260 }; 261 #endif // #if HAVE_ZOLTAN 262 263 enum method_t {HSFC, PHG, PARMETIS}; 264 template< class ldb_vertex_map_t, class ldb_edge_set_t, class ldb_connect_set_t > CALL_Zoltan_LB_Partition(method_t method,MpAccessGlobal & mpa,ldb_vertex_map_t & vertexMap,ldb_edge_set_t & edgeSet,ldb_connect_set_t & connect,const double givenTolerance,const bool verbose)265 bool CALL_Zoltan_LB_Partition( method_t method, 266 MpAccessGlobal &mpa, 267 ldb_vertex_map_t& vertexMap, 268 ldb_edge_set_t& edgeSet, 269 ldb_connect_set_t& connect, 270 const double givenTolerance, 271 const bool verbose ) 272 { 273 #if HAVE_ZOLTAN && HAVE_MPI 274 using ::ALUGrid::MpAccessMPI; 275 276 MpAccessMPI* mpaMPI = dynamic_cast<MpAccessMPI *> (&mpa); 277 if( mpaMPI == 0 ) 278 { 279 std::cerr << "ERROR: wrong mpAccess object, couldn't convert to MpAccessMPI!! in: " << __FILE__ << " line : " << __LINE__ << std::endl; 280 std::abort(); 281 } 282 283 // get communincator (see mpAccess_MPI.cc 284 MPI_Comm comm = mpaMPI->communicator(); 285 int rank = mpaMPI->myrank(); 286 287 typedef ObjectCollection< ldb_vertex_map_t, ldb_edge_set_t > ObjectCollectionType; 288 289 ObjectCollectionType objects( rank, vertexMap, edgeSet ); 290 std::unique_ptr< Zoltan > zz( new Zoltan( comm ) ); 291 alugrid_assert ( zz ); 292 293 // General parameters 294 const char* debug = ( verbose ) ? "1" : "0"; 295 zz->Set_Param( "DEBUG_LEVEL", debug ); 296 zz->Set_Param( "OBJ_WEIGHT_DIM", "1"); 297 zz->Set_Param( "NUM_GID_ENTRIES", "1"); 298 zz->Set_Param( "NUM_LID_ENTRIES", "1"); 299 zz->Set_Param( "RETURN_LISTS", "ALL"); 300 301 if ( method == HSFC ) // edgeSet.size() == 0 ) 302 { 303 zz->Set_Param( "LB_METHOD", "HSFC"); 304 //zz->Set_Param( "KEEP_CUTS", "1" ); 305 zz->Set_Param( "RCB_OUTPUT_LEVEL", "0"); 306 zz->Set_Param( "RCB_RECTILINEAR_BLOCKS", "1"); 307 zz->Set_Param( "LB_APPROACH","REPARTITION"); 308 } 309 else 310 { 311 zz->Set_Param( "LB_METHOD", "GRAPH"); 312 zz->Set_Param( "LB_APPROACH", "REPARTITION"); 313 // zz->Set_Param( "LB_APPROACH","PARTITION"); // give an error in PARMETIS with an 314 // empty partitioning - no idea why... 315 // zz->Set_Param( "LB_APPROACH", "REFINE"); // gives bad loadbalance with 316 // PARMETOS and the rest of the setting - no idea why 317 zz->Set_Param( "EDGE_WEIGHT_DIM","1"); 318 zz->Set_Param( "GRAPH_SYMMETRIZE","NONE" ); 319 // zz->Set_Param( "GRAPH_SYMMETRIZE","TRANSPOSE"); 320 zz->Set_Param( "GRAPH_SYM_WEIGHT","MAX"); 321 // zz->Set_Param( "GRAPH_BUILD_TYPE","FAST_NO_DUP"); 322 #if HAVE_PARMETIS 323 if (method == PARMETIS) 324 zz->Set_Param( "GRAPH_PACKAGE","PARMETIS"); 325 #elif defined(HAVE_SCOTCH) 326 zz->Set_Param( "GRAPH_PACKAGE","SCOTCH"); 327 #endif 328 zz->Set_Param( "CHECK_GRAPH", "0"); 329 zz->Set_Param( "PHG_EDGE_SIZE_THRESHOLD", ".25"); 330 } 331 332 333 334 zz->Set_Num_Obj_Fn ( ObjectCollectionType::get_number_of_objects, &objects); 335 zz->Set_Obj_List_Fn( ObjectCollectionType::get_object_list, &objects); 336 zz->Set_Num_Geom_Fn( ObjectCollectionType::get_num_geometry, &objects); 337 zz->Set_Geom_Multi_Fn( ObjectCollectionType::get_geometry_list, &objects); 338 zz->Set_Num_Edges_Multi_Fn(ObjectCollectionType::get_num_edges_list, &objects); 339 zz->Set_Edge_List_Multi_Fn(ObjectCollectionType::get_edge_list, &objects); 340 341 int changes = 0; 342 int numGidEntries = 0; 343 int numLidEntries = 0; 344 int numImport = 0; 345 ZOLTAN_ID_PTR importGlobalIds = 0; 346 ZOLTAN_ID_PTR importLocalIds = 0; 347 int *importProcs = 0; 348 int *importToPart = 0; 349 int numExport = 0; 350 ZOLTAN_ID_PTR exportGlobalIds = 0; 351 ZOLTAN_ID_PTR exportLocalIds = 0; 352 int *exportProcs = 0; 353 int *exportToPart = 0; 354 355 double tolerance = givenTolerance ; 356 int rc = ZOLTAN_OK + 1 ; 357 int count = 0 ; 358 while( rc != ZOLTAN_OK ) 359 { 360 // std::cout << "[" << rank << "]: "; 361 // std::cout << "zoltan partitioning tolerance: " << tolerance << std::endl; 362 // tolerance for load imbalance 363 { 364 std::stringstream tol; 365 tol << tolerance ; 366 zz->Set_Param( "IMBALANCE_TOL", tol.str() ); 367 } 368 369 rc = zz->LB_Partition(changes, numGidEntries, numLidEntries, 370 numImport, importGlobalIds, importLocalIds, importProcs, importToPart, 371 numExport, exportGlobalIds, exportLocalIds, exportProcs, exportToPart); 372 373 // increase imbalance tolerance and try again 374 if( rc != ZOLTAN_OK ) 375 { 376 tolerance *= 1.05 ; 377 ++ count ; 378 if( count > 3 ) break ; 379 } 380 } 381 382 // if new partitioning has been calculated 383 if (rc == ZOLTAN_OK) 384 { 385 typedef typename ldb_vertex_map_t::iterator iterator; 386 for (int i=0; i < numExport; ++i) 387 { 388 iterator vertex = vertexMap.find( exportGlobalIds[ i ] ); 389 alugrid_assert ( vertex != vertexMap.end () ); 390 (*vertex).second = exportProcs[ i ]; 391 } 392 393 const iterator iEnd = vertexMap.end (); 394 const int myrank = mpaMPI->myrank(); 395 for ( iterator i = vertexMap.begin (); i != iEnd; ++i ) 396 { 397 int& moveTo = i->second; 398 // insert and also set partition number new (including own number) 399 if ( moveTo == -1 ) moveTo = myrank ; 400 connect.insert( MpAccessLocal::sendRank( moveTo ) ); 401 } 402 403 // insert also process number that I will receive objects from 404 for (int i=0; i < numImport; ++i) 405 { 406 connect.insert( MpAccessLocal::recvRank( importProcs[ i ] ) ); 407 } 408 } 409 else 410 { 411 if( verbose && mpa.myrank() == 0 ) 412 std::cerr << "ERROR: Zoltan partitioning failed, partitioning won't change! " << std::endl; 413 // no changes 414 changes = 0; 415 } 416 417 //////////////////////////////////////////////////////////////// 418 // Free the arrays allocated by LB_Partition, and free 419 // the storage allocated for the Zoltan structure and the mesh. 420 //////////////////////////////////////////////////////////////// 421 422 Zoltan::LB_Free_Part(&importGlobalIds, &importLocalIds, &importProcs, &importToPart); 423 Zoltan::LB_Free_Part(&exportGlobalIds, &exportLocalIds, &exportProcs, &exportToPart); 424 425 // delete zoltan structure 426 zz.reset(); 427 428 // return true if partition changed 429 return (changes > 0); 430 #else 431 std::cerr << "ERROR: Zoltan library not found, cannot use Zoltan partitioning! " << std::endl; 432 std::abort(); 433 return false ; 434 #endif // #if HAVE_ZOLTAN && HAVE_MPI 435 } // CALL_Zoltan_LB_Partition 436 437 } // namespace ALUGridZoltan 438 439 #endif // #ifndef ALUGRID_ZOLTAN_H_INCLUDED 440