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