Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
ZoltanPartitioner.cpp
Go to the documentation of this file.
1 /* $Id$
2  *
3  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
4  * storing and accessing finite element mesh data.
5  *
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  *
10  * This library is free software; you can redistribute it and/or
11  * modify it under the terms of the GNU Lesser General Public
12  * License as published by the Free Software Foundation; either
13  * version 2.1 of the License, or (at your option) any later version.
14  *
15  */
16 
17 /** Get a mesh from MOAB and write a Zoltan partition set back into MOAB and to
18  * a file
19  *
20  *
21  */
22 
23 #include <iostream>
24 #include <cassert>
25 #include <sstream>
26 #include <algorithm>
27 
28 #include "moab/ZoltanPartitioner.hpp"
29 #include "moab/Interface.hpp"
30 #include "Internals.hpp"
31 #include "moab/Range.hpp"
32 #include "moab/WriteUtilIface.hpp"
33 #include "moab/MeshTopoUtil.hpp"
34 #include "MBTagConventions.hpp"
35 #include "moab/CN.hpp"
36 // used for gnomonic projection
37 #include "moab/IntxMesh/IntxUtils.hpp"
38 
39 using namespace moab;
40 
41 #define RR \
42  if( MB_SUCCESS != result ) return result
43 
44 static double* Points = NULL;
45 static int* GlobalIds = NULL;
46 static int NumPoints = 0;
47 static int* NumEdges = NULL;
48 static int* NborGlobalId = NULL;
49 static int* NborProcs = NULL;
50 static double* ObjWeights = NULL;
51 static double* EdgeWeights = NULL;
52 static int* Parts = NULL;
53 
54 const bool debug = false;
55 
56 ZoltanPartitioner::ZoltanPartitioner( Interface* impl,
57 #ifdef MOAB_HAVE_MPI
58  ParallelComm* parcomm,
59 #endif
60  const bool use_coords,
61  int argc,
62  char** argv )
63  : PartitionerBase< int >( impl,
64  use_coords
65 #ifdef MOAB_HAVE_MPI
66  ,
67  parcomm
68 #endif
69  ),
70  myZZ( NULL ), myNumPts( 0 ), argcArg( argc ), argvArg( argv )
71 {
72 }
73 
74 ZoltanPartitioner::~ZoltanPartitioner()
75 {
76  if( NULL != myZZ ) delete myZZ;
77 }
78 
79 ErrorCode ZoltanPartitioner::balance_mesh( const char* zmethod,
80  const char* other_method,
81  const bool write_as_sets,
82  const bool write_as_tags )
83 {
84  if( !strcmp( zmethod, "RR" ) && !strcmp( zmethod, "RCB" ) && !strcmp( zmethod, "RIB" ) &&
85  !strcmp( zmethod, "HSFC" ) && !strcmp( zmethod, "Hypergraph" ) && !strcmp( zmethod, "PHG" ) &&
86  !strcmp( zmethod, "PARMETIS" ) && !strcmp( zmethod, "OCTPART" ) )
87  {
88  std::cout << "ERROR node " << mbpc->proc_config().proc_rank() << ": Method must be "
89  << "RR, RCB, RIB, HSFC, Hypergraph (PHG), PARMETIS, or OCTPART" << std::endl;
90  return MB_FAILURE;
91  }
92 
93  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
94  std::vector< int > ids; // point ids from MOAB
95  std::vector< int > adjs, length;
96  Range elems;
97 
98  // Get a mesh from MOAB and divide it across processors.
99 
100  ErrorCode result;
101 
102  if( mbpc->proc_config().proc_rank() == 0 )
103  {
104  result = assemble_graph( 3, pts, ids, adjs, length, elems );RR;
105  }
106 
107  myNumPts = mbInitializePoints( (int)ids.size(), &pts[0], &ids[0], &adjs[0], &length[0] );
108 
109  // Initialize Zoltan. This is a C call. The simple C++ code
110  // that creates Zoltan objects does not keep track of whether
111  // Zoltan_Initialize has been called.
112 
113  float version;
114 
115  Zoltan_Initialize( argcArg, argvArg, &version );
116 
117  // Create Zoltan object. This calls Zoltan_Create.
118  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
119 
120  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
121  SetRCB_Parameters();
122  else if( !strcmp( zmethod, "RIB" ) )
123  SetRIB_Parameters();
124  else if( !strcmp( zmethod, "HSFC" ) )
125  SetHSFC_Parameters();
126  else if( !strcmp( zmethod, "Hypergraph" ) || !strcmp( zmethod, "PHG" ) )
127  if( NULL == other_method )
128  SetHypergraph_Parameters( "auto" );
129  else
130  SetHypergraph_Parameters( other_method );
131  else if( !strcmp( zmethod, "PARMETIS" ) )
132  {
133  if( NULL == other_method )
134  SetPARMETIS_Parameters( "RepartGDiffusion" );
135  else
136  SetPARMETIS_Parameters( other_method );
137  }
138  else if( !strcmp( zmethod, "OCTPART" ) )
139  {
140  if( NULL == other_method )
141  SetOCTPART_Parameters( "2" );
142  else
143  SetOCTPART_Parameters( other_method );
144  }
145 
146  // Call backs:
147 
148  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
149  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
150  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
151  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
152  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
153  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
154 
155  // Perform the load balancing partitioning
156 
157  int changes;
158  int numGidEntries;
159  int numLidEntries;
160  int numImport;
161  ZOLTAN_ID_PTR importGlobalIds;
162  ZOLTAN_ID_PTR importLocalIds;
163  int* importProcs;
164  int* importToPart;
165  int numExport;
166  ZOLTAN_ID_PTR exportGlobalIds;
167  ZOLTAN_ID_PTR exportLocalIds;
168  int* exportProcs;
169  int* exportToPart;
170 
171  int rc = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, numImport, importGlobalIds, importLocalIds,
172  importProcs, importToPart, numExport, exportGlobalIds, exportLocalIds, exportProcs,
173  exportToPart );
174 
175  rc = mbGlobalSuccess( rc );
176 
177  if( !rc )
178  {
179  mbPrintGlobalResult( "==============Result==============", myNumPts, numImport, numExport, changes );
180  }
181  else
182  {
183  return MB_FAILURE;
184  }
185 
186  // take results & write onto MOAB partition sets
187 
188  int* assignment;
189 
190  mbFinalizePoints( (int)ids.size(), numExport, exportLocalIds, exportProcs, &assignment );
191 
192  if( mbpc->proc_config().proc_rank() == 0 )
193  {
194  result = write_partition( mbpc->proc_config().proc_size(), elems, assignment, write_as_sets, write_as_tags );
195 
196  if( MB_SUCCESS != result ) return result;
197 
198  free( (int*)assignment );
199  }
200 
201  // Free the memory allocated for lists returned by LB_Parition()
202 
203  myZZ->LB_Free_Part( &importGlobalIds, &importLocalIds, &importProcs, &importToPart );
204  myZZ->LB_Free_Part( &exportGlobalIds, &exportLocalIds, &exportProcs, &exportToPart );
205 
206  // Implementation note: A Zoltan object contains an MPI communicator.
207  // When the Zoltan object is destroyed, it uses it's MPI communicator.
208  // So it is important that the Zoltan object is destroyed before
209  // the MPI communicator is destroyed. To ensure this, dynamically
210  // allocate the Zoltan object, so you can explicitly destroy it.
211  // If you create a Zoltan object on the stack, it's destructor will
212  // be invoked atexit, possibly after the communicator's
213  // destructor.
214 
215  return MB_SUCCESS;
216 }
217 ErrorCode ZoltanPartitioner::repartition_to_procs( std::vector< double >& x,
218  std::vector< double >& y,
219  std::vector< double >& z,
220  std::vector< int >& ids,
221  const char* zmethod,
222  std::vector< int >& dest )
223 {
224  //
225  int nprocs = mbpc->proc_config().proc_size();
226  int rank = mbpc->proc_config().proc_rank();
227  clock_t t = clock();
228  // form pts and ids, as in assemble_graph
229  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
230  pts.resize( x.size() * 3 );
231  for( size_t i = 0; i < x.size(); i++ )
232  {
233  pts[3 * i] = x[i];
234  pts[3 * i + 1] = y[i];
235  pts[3 * i + 2] = z[i];
236  }
237  // do not get out until done!
238 
239  Points = &pts[0];
240  GlobalIds = &ids[0];
241  NumPoints = (int)x.size();
242  NumEdges = NULL;
243  NborGlobalId = NULL;
244  NborProcs = NULL;
245  ObjWeights = NULL;
246  EdgeWeights = NULL;
247  Parts = NULL;
248 
249  float version;
250  if( rank == 0 ) std::cout << "Initializing zoltan..." << std::endl;
251 
252  Zoltan_Initialize( argcArg, argvArg, &version );
253 
254  // Create Zoltan object. This calls Zoltan_Create.
255  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
256 
257  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
258  SetRCB_Parameters();
259  else if( !strcmp( zmethod, "RIB" ) )
260  SetRIB_Parameters();
261  else if( !strcmp( zmethod, "HSFC" ) )
262  SetHSFC_Parameters();
263 
264  // set # requested partitions
265  char buff[10];
266  snprintf( buff, 10, "%d", nprocs );
267  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
268  if( ZOLTAN_OK != retval ) return MB_FAILURE;
269 
270  // request all, import and export
271  retval = myZZ->Set_Param( "RETURN_LISTS", "ALL" );
272  if( ZOLTAN_OK != retval ) return MB_FAILURE;
273 
274  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
275  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
276  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
277  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
278  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
279  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
280 
281  // Perform the load balancing partitioning
282 
283  int changes;
284  int numGidEntries;
285  int numLidEntries;
286  int num_import;
287  ZOLTAN_ID_PTR import_global_ids, import_local_ids;
288  int* import_procs;
289  int* import_to_part;
290  int num_export;
291  ZOLTAN_ID_PTR export_global_ids, export_local_ids;
292  int *assign_procs, *assign_parts;
293 
294  if( rank == 0 )
295  std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nprocs
296  << " processors..." << std::endl;
297 
298  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids,
299  import_procs, import_to_part, num_export, export_global_ids, export_local_ids,
300  assign_procs, assign_parts );
301 
302  if( ZOLTAN_OK != retval ) return MB_FAILURE;
303 
304  if( rank == 0 )
305  {
306  std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
307  t = clock();
308  }
309  for( int i = 0; i < NumPoints; i++ )
310  dest[i] = rank;
311  for( int i = 0; i < num_export; i++ )
312  dest[export_local_ids[i]] = assign_procs[i];
313 
314  // Free data structures allocated by Zoltan::LB_Partition
315  retval = myZZ->LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part );
316  if( ZOLTAN_OK != retval ) return MB_FAILURE;
317  retval = myZZ->LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts );
318  if( ZOLTAN_OK != retval ) return MB_FAILURE;
319  return MB_SUCCESS;
320 }
321 
322 ErrorCode ZoltanPartitioner::repartition( std::vector< double >& x,
323  std::vector< double >& y,
324  std::vector< double >& z,
325  int StartID,
326  const char* zmethod,
327  Range& localGIDs )
328 {
329  //
330  int nprocs = mbpc->proc_config().proc_size();
331  int rank = mbpc->proc_config().proc_rank();
332  clock_t t = clock();
333  // form pts and ids, as in assemble_graph
334  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
335  pts.resize( x.size() * 3 );
336  std::vector< int > ids; // point ids from MOAB
337  ids.resize( x.size() );
338  for( size_t i = 0; i < x.size(); i++ )
339  {
340  pts[3 * i] = x[i];
341  pts[3 * i + 1] = y[i];
342  pts[3 * i + 2] = z[i];
343  ids[i] = StartID + (int)i;
344  }
345  // do not get out until done!
346 
347  Points = &pts[0];
348  GlobalIds = &ids[0];
349  NumPoints = (int)x.size();
350  NumEdges = NULL;
351  NborGlobalId = NULL;
352  NborProcs = NULL;
353  ObjWeights = NULL;
354  EdgeWeights = NULL;
355  Parts = NULL;
356 
357  float version;
358  if( rank == 0 ) std::cout << "Initializing zoltan..." << std::endl;
359 
360  Zoltan_Initialize( argcArg, argvArg, &version );
361 
362  // Create Zoltan object. This calls Zoltan_Create.
363  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
364 
365  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
366  SetRCB_Parameters();
367  else if( !strcmp( zmethod, "RIB" ) )
368  SetRIB_Parameters();
369  else if( !strcmp( zmethod, "HSFC" ) )
370  SetHSFC_Parameters();
371 
372  // set # requested partitions
373  char buff[10];
374  snprintf( buff, 10, "%d", nprocs );
375  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
376  if( ZOLTAN_OK != retval ) return MB_FAILURE;
377 
378  // request all, import and export
379  retval = myZZ->Set_Param( "RETURN_LISTS", "ALL" );
380  if( ZOLTAN_OK != retval ) return MB_FAILURE;
381 
382  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
383  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
384  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
385  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
386  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
387  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
388 
389  // Perform the load balancing partitioning
390 
391  int changes;
392  int numGidEntries;
393  int numLidEntries;
394  int num_import;
395  ZOLTAN_ID_PTR import_global_ids, import_local_ids;
396  int* import_procs;
397  int* import_to_part;
398  int num_export;
399  ZOLTAN_ID_PTR export_global_ids, export_local_ids;
400  int *assign_procs, *assign_parts;
401 
402  if( rank == 0 )
403  std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nprocs
404  << " processors..." << std::endl;
405 
406  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids,
407  import_procs, import_to_part, num_export, export_global_ids, export_local_ids,
408  assign_procs, assign_parts );
409  if( ZOLTAN_OK != retval ) return MB_FAILURE;
410 
411  if( rank == 0 )
412  {
413  std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
414  t = clock();
415  }
416 
417  std::sort( import_global_ids, import_global_ids + num_import, std::greater< int >() );
418  std::sort( export_global_ids, export_global_ids + num_export, std::greater< int >() );
419 
420  Range iniGids( (EntityHandle)StartID, (EntityHandle)StartID + x.size() - 1 );
421  Range imported, exported;
422  std::copy( import_global_ids, import_global_ids + num_import, range_inserter( imported ) );
423  std::copy( export_global_ids, export_global_ids + num_export, range_inserter( exported ) );
424  localGIDs = subtract( iniGids, exported );
425  localGIDs = unite( localGIDs, imported );
426  // Free data structures allocated by Zoltan::LB_Partition
427  retval = myZZ->LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part );
428  if( ZOLTAN_OK != retval ) return MB_FAILURE;
429  retval = myZZ->LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts );
430  if( ZOLTAN_OK != retval ) return MB_FAILURE;
431 
432  return MB_SUCCESS;
433 }
434 
435 ErrorCode ZoltanPartitioner::partition_inferred_mesh( EntityHandle sfileset,
436  size_t num_parts,
437  int part_dim,
438  const bool write_as_sets,
439  int projection_type )
440 {
441  ErrorCode result;
442 
443  moab::Range elverts;
444  result = mbImpl->get_entities_by_dimension( sfileset, part_dim, elverts );RR;
445 
446  std::vector< double > elcoords( elverts.size() * 3 );
447  result = mbImpl->get_coords( elverts, &elcoords[0] );RR;
448 
449  std::vector< std::vector< EntityHandle > > part_assignments( num_parts );
450  int part, proc;
451 
452  // Loop over coordinates
453  for( size_t iel = 0; iel < elverts.size(); iel++ )
454  {
455  // Gather coordinates into temporary array
456  double* ecoords = &elcoords[iel * 3];
457 
458  // do a projection if needed
459  if( projection_type > 0 ) IntxUtils::transform_coordinates( ecoords, projection_type );
460 
461  // Compute the coordinate's part assignment
462  myZZ->LB_Point_PP_Assign( ecoords, proc, part );
463 
464  // Store the part assignment in the return array
465  part_assignments[part].push_back( elverts[iel] );
466  }
467 
468  // get the partition set tag
469  Tag part_set_tag = mbpc->partition_tag();
470 
471  // If some entity sets exist, let us delete those first.
472  if( write_as_sets )
473  {
474  moab::Range oldpsets;
475  result = mbImpl->get_entities_by_type_and_tag( sfileset, MBENTITYSET, &part_set_tag, NULL, 1, oldpsets,
476  Interface::UNION );RR;
477  result = mbImpl->remove_entities( sfileset, oldpsets );RR;
478  }
479 
480  size_t nparts_assigned = 0;
481  for( size_t partnum = 0; partnum < num_parts; ++partnum )
482  {
483  std::vector< EntityHandle >& partvec = part_assignments[partnum];
484 
485  nparts_assigned += ( partvec.size() ? 1 : 0 );
486 
487  if( write_as_sets )
488  {
489  EntityHandle partNset;
490  result = mbImpl->create_meshset( moab::MESHSET_SET, partNset );RR;
491  if( partvec.size() )
492  {
493  result = mbImpl->add_entities( partNset, &partvec[0], partvec.size() );RR;
494  }
495  result = mbImpl->add_parent_child( sfileset, partNset );RR;
496 
497  // part numbering should start from 0
498  int ipartnum = (int)partnum;
499 
500  // assign partitions as a sparse tag by grouping elements under sets
501  result = mbImpl->tag_set_data( part_set_tag, &partNset, 1, &ipartnum );RR;
502  }
503  else
504  {
505  /* assign as a dense vector to all elements
506  allocate integer-size partitions */
507  std::vector< int > assignment( partvec.size(),
508  partnum ); // initialize all values to partnum
509  result = mbImpl->tag_set_data( part_set_tag, &partvec[0], partvec.size(), &assignment[0] );RR;
510  }
511  }
512 
513  if( nparts_assigned != num_parts )
514  {
515  std::cout << "WARNING: The inference yielded lesser number of parts (" << nparts_assigned
516  << ") than requested by user (" << num_parts << ").\n";
517  }
518 
519  return MB_SUCCESS;
520 }
521 
522 ErrorCode ZoltanPartitioner::partition_mesh_and_geometry( const double part_geom_mesh_size,
523  const int nparts,
524  const char* zmethod,
525  const char* other_method,
526  double imbal_tol,
527  const int part_dim,
528  const bool write_as_sets,
529  const bool write_as_tags,
530  const int obj_weight,
531  const int edge_weight,
532  const int projection_type,
533  const bool recompute_rcb_box,
534  const bool print_time )
535 {
536  // should only be called in serial
537  if( mbpc->proc_config().proc_size() != 1 )
538  {
539  std::cout << "ZoltanPartitioner::partition_mesh_and_geometry must be called in serial." << std::endl;
540  return MB_FAILURE;
541  }
542  clock_t t = clock();
543  if( NULL != zmethod && strcmp( zmethod, "RR" ) && strcmp( zmethod, "RCB" ) && strcmp( zmethod, "RIB" ) &&
544  strcmp( zmethod, "HSFC" ) && strcmp( zmethod, "Hypergraph" ) && strcmp( zmethod, "PHG" ) &&
545  strcmp( zmethod, "PARMETIS" ) && strcmp( zmethod, "OCTPART" ) )
546  {
547  std::cout << "ERROR node " << mbpc->proc_config().proc_rank() << ": Method must be "
548  << "RCB, RIB, HSFC, Hypergraph (PHG), PARMETIS, or OCTPART" << std::endl;
549  return MB_FAILURE;
550  }
551 
552  bool part_geom = false;
553  if( 0 == strcmp( zmethod, "RR" ) || 0 == strcmp( zmethod, "RCB" ) || 0 == strcmp( zmethod, "RIB" ) ||
554  0 == strcmp( zmethod, "HSFC" ) )
555  part_geom = true; // so no adjacency / edges needed
556  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
557  std::vector< int > ids; // point ids from MOAB
558  std::vector< int > adjs, length, parts;
559  std::vector< double > obj_weights, edge_weights;
560  Range elems;
561 
562  // Get a mesh from MOAB and diide it across processors.
563 
564  ErrorCode result = MB_SUCCESS;
565 
566  // short-circuit everything if RR partition is requested
567  if( !strcmp( zmethod, "RR" ) )
568  {
569  if( part_geom_mesh_size < 0. )
570  {
571  // get all elements
572  result = mbImpl->get_entities_by_dimension( 0, part_dim, elems );RR;
573  if( elems.empty() ) return MB_FAILURE;
574  // make a trivial assignment vector
575  std::vector< int > assign_vec( elems.size() );
576  int num_per = elems.size() / nparts;
577  int extra = elems.size() % nparts;
578  if( extra ) num_per++;
579  int nstart = 0;
580  for( int i = 0; i < nparts; i++ )
581  {
582  if( i == extra ) num_per--;
583  std::fill( &assign_vec[nstart], &assign_vec[nstart + num_per], i );
584  nstart += num_per;
585  }
586 
587  result = write_partition( nparts, elems, &assign_vec[0], write_as_sets, write_as_tags );RR;
588  }
589  else
590  {
591  }
592 
593  return result;
594  }
595 
596  std::cout << "Assembling graph..." << std::endl;
597 
598  if( part_geom_mesh_size < 0. )
599  {
600  // if (!part_geom) {
601  result = assemble_graph( part_dim, pts, ids, adjs, length, elems, part_geom, projection_type );RR;
602  }
603  else
604  {
605  MB_CHK_SET_ERR( MB_FAILURE, "Geometry partitions not supported.\n" );
606  }
607  if( print_time )
608  {
609  std::cout << " time to assemble graph: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
610  t = clock();
611  }
612  double* o_wgt = NULL;
613  double* e_wgt = NULL;
614  if( obj_weights.size() > 0 ) o_wgt = &obj_weights[0];
615  if( edge_weights.size() > 0 ) e_wgt = &edge_weights[0];
616 
617  myNumPts = mbInitializePoints( (int)ids.size(), &pts[0], &ids[0], &adjs[0], &length[0], o_wgt, e_wgt, &parts[0],
618  part_geom );
619 
620  // Initialize Zoltan. This is a C call. The simple C++ code
621  // that creates Zoltan objects does not keep track of whether
622  // Zoltan_Initialize has been called.
623 
624  if( print_time )
625  {
626  std::cout << " time to initialize points: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
627  t = clock();
628  }
629  float version;
630 
631  std::cout << "Initializing zoltan..." << std::endl;
632 
633  Zoltan_Initialize( argcArg, argvArg, &version );
634 
635  // Create Zoltan object. This calls Zoltan_Create.
636  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
637 
638  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
639  SetRCB_Parameters( recompute_rcb_box );
640  else if( !strcmp( zmethod, "RIB" ) )
641  SetRIB_Parameters();
642  else if( !strcmp( zmethod, "HSFC" ) )
643  SetHSFC_Parameters();
644  else if( !strcmp( zmethod, "Hypergraph" ) || !strcmp( zmethod, "PHG" ) )
645  {
646  if( NULL == other_method || ( other_method[0] == '\0' ) )
647  SetHypergraph_Parameters( "auto" );
648  else
649  SetHypergraph_Parameters( other_method );
650 
651  if( imbal_tol )
652  {
653  std::ostringstream str;
654  str << imbal_tol;
655  myZZ->Set_Param( "IMBALANCE_TOL", str.str().c_str() ); // no debug messages
656  }
657  }
658  else if( !strcmp( zmethod, "PARMETIS" ) )
659  {
660  if( NULL == other_method )
661  SetPARMETIS_Parameters( "RepartGDiffusion" );
662  else
663  SetPARMETIS_Parameters( other_method );
664  }
665  else if( !strcmp( zmethod, "OCTPART" ) )
666  {
667  if( NULL == other_method )
668  SetOCTPART_Parameters( "2" );
669  else
670  SetOCTPART_Parameters( other_method );
671  }
672 
673  // set # requested partitions
674  char buff[10];
675  snprintf( buff, 10, "%d", nparts );
676  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
677  if( ZOLTAN_OK != retval ) return MB_FAILURE;
678 
679  // request only partition assignments
680  retval = myZZ->Set_Param( "RETURN_LISTS", "PARTITION ASSIGNMENTS" );
681  if( ZOLTAN_OK != retval ) return MB_FAILURE;
682 
683  if( obj_weight > 0 )
684  {
685  std::ostringstream str;
686  str << obj_weight;
687  retval = myZZ->Set_Param( "OBJ_WEIGHT_DIM", str.str().c_str() );
688  if( ZOLTAN_OK != retval ) return MB_FAILURE;
689  }
690 
691  if( edge_weight > 0 )
692  {
693  std::ostringstream str;
694  str << edge_weight;
695  retval = myZZ->Set_Param( "EDGE_WEIGHT_DIM", str.str().c_str() );
696  if( ZOLTAN_OK != retval ) return MB_FAILURE;
697  }
698 
699  // Call backs:
700 
701  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
702  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
703  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
704  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
705  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
706  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
707  if( part_geom_mesh_size > 0. )
708  {
709  myZZ->Set_Part_Multi_Fn( mbGetPart, NULL );
710  }
711 
712  // Perform the load balancing partitioning
713 
714  int changes;
715  int numGidEntries;
716  int numLidEntries;
717  int dumnum1;
718  ZOLTAN_ID_PTR dum_local, dum_global;
719  int *dum1, *dum2;
720  int num_assign;
721  ZOLTAN_ID_PTR assign_gid, assign_lid;
722  int *assign_procs, *assign_parts;
723 
724  std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nparts
725  << " processors..." << std::endl;
726 #ifndef NDEBUG
727 #if 0
728  if (NULL == zmethod || !strcmp(zmethod, "RCB"))
729  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)zmethod, 1, 1, 0, 0);
730 
731  if ( !strcmp(zmethod, "PHG"))
732  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)zmethod, 1, 0, 1, 1);
733 #endif
734 #endif
735  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, dumnum1, dum_global, dum_local, dum1, dum2,
736  num_assign, assign_gid, assign_lid, assign_procs, assign_parts );
737  if( ZOLTAN_OK != retval ) return MB_FAILURE;
738 
739  if( print_time )
740  {
741  std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
742  t = clock();
743  }
744  // take results & write onto MOAB partition sets
745  std::cout << "Saving partition information to MOAB..." << std::endl;
746 
747  if( part_geom_mesh_size < 0. )
748  {
749  // if (!part_geom) {
750  result = write_partition( nparts, elems, assign_parts, write_as_sets, write_as_tags );
751  }
752  else
753  {
754  MB_CHK_SET_ERR( MB_FAILURE, "Geometry partitions not supported.\n" );
755  }
756 
757  if( print_time )
758  {
759  std::cout << " time to write partition in memory " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
760  t = clock();
761  }
762 
763  if( MB_SUCCESS != result ) return result;
764 
765  // Free the memory allocated for lists returned by LB_Parition()
766  myZZ->LB_Free_Part( &assign_gid, &assign_lid, &assign_procs, &assign_parts );
767 
768  return MB_SUCCESS;
769 }
770 
771 ErrorCode ZoltanPartitioner::include_closure()
772 {
773  ErrorCode result;
774  Range ents;
775  Range adjs;
776  std::cout << "Adding closure..." << std::endl;
777 
778  for( Range::iterator rit = partSets.begin(); rit != partSets.end(); ++rit )
779  {
780 
781  // get the top-dimensional entities in the part
782  result = mbImpl->get_entities_by_handle( *rit, ents, true );RR;
783 
784  if( ents.empty() ) continue;
785 
786  // get intermediate-dimensional adjs and add to set
787  for( int d = mbImpl->dimension_from_handle( *ents.begin() ) - 1; d >= 1; d-- )
788  {
789  adjs.clear();
790  result = mbImpl->get_adjacencies( ents, d, false, adjs, Interface::UNION );RR;
791  result = mbImpl->add_entities( *rit, adjs );RR;
792  }
793 
794  // now get vertices and add to set; only need to do for ents, not for adjs
795  adjs.clear();
796  result = mbImpl->get_adjacencies( ents, 0, false, adjs, Interface::UNION );RR;
797  result = mbImpl->add_entities( *rit, adjs );RR;
798 
799  ents.clear();
800  }
801 
802  // now go over non-part entity sets, looking for contained entities
803  Range sets, part_ents;
804  result = mbImpl->get_entities_by_type( 0, MBENTITYSET, sets );RR;
805  for( Range::iterator rit = sets.begin(); rit != sets.end(); ++rit )
806  {
807  // skip parts
808  if( partSets.find( *rit ) != partSets.end() ) continue;
809 
810  // get entities in this set, recursively
811  ents.clear();
812  result = mbImpl->get_entities_by_handle( *rit, ents, true );RR;
813 
814  // now check over all parts
815  for( Range::iterator rit2 = partSets.begin(); rit2 != partSets.end(); ++rit2 )
816  {
817  part_ents.clear();
818  result = mbImpl->get_entities_by_handle( *rit2, part_ents, false );RR;
819  Range int_range = intersect( ents, part_ents );
820  if( !int_range.empty() )
821  {
822  // non-empty intersection, add to part set
823  result = mbImpl->add_entities( *rit2, &( *rit ), 1 );RR;
824  }
825  }
826  }
827 
828  // finally, mark all the part sets as having the closure
829  Tag closure_tag;
830  result =
831  mbImpl->tag_get_handle( "INCLUDES_CLOSURE", 1, MB_TYPE_INTEGER, closure_tag, MB_TAG_SPARSE | MB_TAG_CREAT );RR;
832 
833  std::vector< int > closure_vals( partSets.size(), 1 );
834  result = mbImpl->tag_set_data( closure_tag, partSets, &closure_vals[0] );RR;
835 
836  return MB_SUCCESS;
837 }
838 
839 ErrorCode ZoltanPartitioner::assemble_graph( const int dimension,
840  std::vector< double >& coords,
841  std::vector< int >& moab_ids,
842  std::vector< int >& adjacencies,
843  std::vector< int >& length,
844  Range& elems,
845  bool part_geom,
846  int projection_type )
847 {
848  // assemble a graph with vertices equal to elements of specified dimension, edges
849  // signified by list of other elements to which an element is connected
850 
851  // get the elements of that dimension
852  ErrorCode result = mbImpl->get_entities_by_dimension( 0, dimension, elems );
853  if( MB_SUCCESS != result || elems.empty() ) return result;
854 
855  // assign global ids
856  if( assign_global_ids )
857  {
858  EntityHandle rootset = 0;
859  result = mbpc->assign_global_ids( rootset, dimension, 1, true, true );RR;
860  }
861 
862  // now assemble the graph, calling MeshTopoUtil to get bridge adjacencies through d-1
863  // dimensional neighbors
864  MeshTopoUtil mtu( mbImpl );
865  Range adjs;
866  // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited
867  // by MBCN
868  int neighbors[5 * MAX_SUB_ENTITIES];
869  double avg_position[3];
870  int moab_id;
871 
872  // get the global id tag handle
873  Tag gid = mbImpl->globalId_tag();
874 
875  for( Range::iterator rit = elems.begin(); rit != elems.end(); ++rit )
876  {
877 
878  if( !part_geom )
879  {
880  // get bridge adjacencies
881  adjs.clear();
882  result = mtu.get_bridge_adjacencies( *rit, ( dimension > 0 ? dimension - 1 : 3 ), dimension, adjs );RR;
883 
884  // get the graph vertex ids of those
885  if( !adjs.empty() )
886  {
887  assert( adjs.size() < 5 * MAX_SUB_ENTITIES );
888  result = mbImpl->tag_get_data( gid, adjs, neighbors );RR;
889  }
890 
891  // copy those into adjacencies vector
892  length.push_back( (int)adjs.size() );
893  std::copy( neighbors, neighbors + adjs.size(), std::back_inserter( adjacencies ) );
894  }
895 
896  // get average position of vertices
897  result = mtu.get_average_position( *rit, avg_position );RR;
898 
899  // get the graph vertex id for this element
900  result = mbImpl->tag_get_data( gid, &( *rit ), 1, &moab_id );RR;
901 
902  // copy those into coords vector
903  moab_ids.push_back( moab_id );
904  // transform coordinates to spherical coordinates, if requested
905  if( projection_type > 0 ) IntxUtils::transform_coordinates( avg_position, projection_type );
906 
907  std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) );
908  }
909 
910  if( debug )
911  {
912  std::cout << "Length vector: " << std::endl;
913  std::copy( length.begin(), length.end(), std::ostream_iterator< int >( std::cout, ", " ) );
914  std::cout << std::endl;
915  std::cout << "Adjacencies vector: " << std::endl;
916  std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< int >( std::cout, ", " ) );
917  std::cout << std::endl;
918  std::cout << "Moab_ids vector: " << std::endl;
919  std::copy( moab_ids.begin(), moab_ids.end(), std::ostream_iterator< int >( std::cout, ", " ) );
920  std::cout << std::endl;
921  std::cout << "Coords vector: " << std::endl;
922  std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( std::cout, ", " ) );
923  std::cout << std::endl;
924  }
925 
926  return MB_SUCCESS;
927 }
928 
929 ErrorCode ZoltanPartitioner::write_partition( const int nparts,
930  Range& elems,
931  const int* assignment,
932  const bool write_as_sets,
933  const bool write_as_tags )
934 {
935  ErrorCode result;
936 
937  // get the partition set tag
938  Tag part_set_tag;
939  int dum_id = -1, i;
940  result = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag,
941  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );RR;
942 
943  // get any sets already with this tag, and clear them
944  Range tagged_sets;
945  result =
946  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );RR;
947  if( !tagged_sets.empty() )
948  {
949  result = mbImpl->clear_meshset( tagged_sets );RR;
950  if( !write_as_sets )
951  {
952  result = mbImpl->tag_delete_data( part_set_tag, tagged_sets );RR;
953  }
954  }
955 
956  if( write_as_sets )
957  {
958  // first, create partition sets and store in vector
959  partSets.clear();
960 
961  if( nparts > (int)tagged_sets.size() )
962  {
963  // too few partition sets - create missing ones
964  int num_new = nparts - tagged_sets.size();
965  for( i = 0; i < num_new; i++ )
966  {
967  EntityHandle new_set;
968  result = mbImpl->create_meshset( MESHSET_SET, new_set );RR;
969  tagged_sets.insert( new_set );
970  }
971  }
972  else if( nparts < (int)tagged_sets.size() )
973  {
974  // too many partition sets - delete extras
975  int num_del = tagged_sets.size() - nparts;
976  for( i = 0; i < num_del; i++ )
977  {
978  EntityHandle old_set = tagged_sets.pop_back();
979  result = mbImpl->delete_entities( &old_set, 1 );RR;
980  }
981  }
982  // assign partition sets to vector
983  partSets.swap( tagged_sets );
984 
985  // write a tag to those sets denoting they're partition sets, with a value of the
986  // proc number
987  int* dum_ids = new int[nparts];
988  for( i = 0; i < nparts; i++ )
989  dum_ids[i] = i;
990 
991  result = mbImpl->tag_set_data( part_set_tag, partSets, dum_ids );RR;
992  // found out by valgrind when we run mbpart
993  delete[] dum_ids;
994  dum_ids = NULL;
995 
996  // assign entities to the relevant sets
997  std::vector< EntityHandle > tmp_part_sets;
998  // int N = (int)elems.size();
999  std::copy( partSets.begin(), partSets.end(), std::back_inserter( tmp_part_sets ) );
1000  /*Range::reverse_iterator riter;
1001  for (i = N-1, riter = elems.rbegin(); riter != elems.rend(); ++riter, i--) {
1002  int assigned_part = assignment[i];
1003  part_ranges[assigned_part].insert(*riter);
1004  //result = mbImpl->add_entities(tmp_part_sets[assignment[i]], &(*rit), 1); RR;
1005  }*/
1006 
1007  Range::iterator rit;
1008  for( i = 0, rit = elems.begin(); rit != elems.end(); ++rit, i++ )
1009  {
1010  result = mbImpl->add_entities( tmp_part_sets[assignment[i]], &( *rit ), 1 );RR;
1011  }
1012  /*for (i=0; i<nparts; i++)
1013  {
1014  result = mbImpl->add_entities(tmp_part_sets[i], part_ranges[i]); RR;
1015  }
1016  delete [] part_ranges;*/
1017  // check for empty sets, warn if there are any
1018  Range empty_sets;
1019  for( rit = partSets.begin(); rit != partSets.end(); ++rit )
1020  {
1021  int num_ents = 0;
1022  result = mbImpl->get_number_entities_by_handle( *rit, num_ents );
1023  if( MB_SUCCESS != result || !num_ents ) empty_sets.insert( *rit );
1024  }
1025  if( !empty_sets.empty() )
1026  {
1027  std::cout << "WARNING: " << empty_sets.size() << " empty sets in partition: ";
1028  for( rit = empty_sets.begin(); rit != empty_sets.end(); ++rit )
1029  std::cout << *rit << " ";
1030  std::cout << std::endl;
1031  }
1032  }
1033 
1034  if( write_as_tags )
1035  {
1036  // allocate integer-size partitions
1037  result = mbImpl->tag_set_data( part_set_tag, elems, assignment );RR;
1038  }
1039 
1040  return MB_SUCCESS;
1041 }
1042 
1043 void ZoltanPartitioner::SetRCB_Parameters( const bool recompute_rcb_box )
1044 {
1045  if( mbpc->proc_config().proc_rank() == 0 )
1046  std::cout << "Using Zoltan-RCB (Recursive Coordinate Bisection) repartitioning scheme..." << std::endl;
1047 
1048  // General parameters:
1049  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1050  myZZ->Set_Param( "LB_METHOD", "RCB" ); // recursive coordinate bisection
1051 
1052  // RCB parameters:
1053  myZZ->Set_Param( "RCB_OUTPUT_LEVEL", "0" ); // increase to 1 for verbose details
1054  myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition so that we can infer partitions
1055  // myZZ->Set_Param("RCB_RECTILINEAR_BLOCKS", "1"); // don't split point on boundary
1056  if( recompute_rcb_box ) myZZ->Set_Param( "RCB_RECOMPUTE_BOX", "1" );
1057 }
1058 
1059 void ZoltanPartitioner::SetRIB_Parameters()
1060 {
1061  if( mbpc->proc_config().proc_rank() == 0 )
1062  std::cout << "Using Zoltan-RIB (Recursive Inertial Bisection) repartitioning scheme..." << std::endl;
1063 
1064  // General parameters:
1065  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1066  myZZ->Set_Param( "LB_METHOD", "RIB" ); // Recursive Inertial Bisection
1067 
1068  // RIB parameters:
1069  myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition
1070  myZZ->Set_Param( "AVERAGE_CUTS", "1" );
1071 }
1072 
1073 void ZoltanPartitioner::SetHSFC_Parameters()
1074 {
1075  if( mbpc->proc_config().proc_rank() == 0 )
1076  std::cout << "Using Zoltan-HSFC (Hilbert Space Filling Curve) repartitioning scheme..." << std::endl;
1077 
1078  // General parameters:
1079  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1080  myZZ->Set_Param( "LB_METHOD", "HSFC" ); // perform Hilbert space filling curve
1081 
1082  // HSFC parameters:
1083  myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition
1084 }
1085 
1086 void ZoltanPartitioner::SetHypergraph_Parameters( const char* phg_method )
1087 {
1088  if( mbpc->proc_config().proc_rank() == 0 )
1089  std::cout << "Using Zoltan-PHG (Hypergraph) repartitioning scheme..." << std::endl;
1090 
1091  // General parameters:
1092  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1093  myZZ->Set_Param( "LB_METHOD", "Hypergraph" ); // Hypergraph (or PHG)
1094 
1095  // Hypergraph (or PHG) parameters:
1096  myZZ->Set_Param( "PHG_COARSEPARTITION_METHOD", phg_method ); // CoarsePartitionMethod
1097 }
1098 
1099 void ZoltanPartitioner::SetPARMETIS_Parameters( const char* parmetis_method )
1100 {
1101  if( mbpc->proc_config().proc_rank() == 0 )
1102  std::cout << "Using Zoltan-ParMetis (" << parmetis_method << ") repartitioning scheme..." << std::endl;
1103 
1104  // General parameters:
1105  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1106  myZZ->Set_Param( "LB_METHOD", "PARMETIS" ); // the ParMETIS library
1107 
1108  // PARMETIS parameters:
1109  myZZ->Set_Param( "PARMETIS_METHOD", parmetis_method ); // method in the library
1110 }
1111 
1112 void ZoltanPartitioner::SetOCTPART_Parameters( const char* oct_method )
1113 {
1114  if( mbpc->proc_config().proc_rank() == 0 )
1115  std::cout << "Using Zoltan-OctTree (" << oct_method << ") repartitioning scheme..." << std::endl;
1116 
1117  // General parameters:
1118  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1119  myZZ->Set_Param( "LB_METHOD", "OCTPART" ); // octree partitioning
1120 
1121  // OCTPART parameters:
1122  myZZ->Set_Param( "OCT_METHOD", oct_method ); // the SFC to be used
1123  myZZ->Set_Param( "OCT_OUTPUT_LEVEL", "3" );
1124 }
1125 
1126 int ZoltanPartitioner::mbInitializePoints( int npts,
1127  double* pts,
1128  int* ids,
1129  int* adjs,
1130  int* length,
1131  double* obj_weights,
1132  double* edge_weights,
1133  int* parts,
1134  bool part_geom )
1135 {
1136  unsigned int i;
1137  int j;
1138  int *numPts, *nborProcs = NULL;
1139  int sum, ptsPerProc, ptsAssigned, mySize;
1140  MPI_Status stat;
1141  double* sendPts;
1142  int* sendIds;
1143  int* sendEdges = NULL;
1144  int* sendNborId = NULL;
1145  int* sendProcs;
1146 
1147  if( mbpc->proc_config().proc_rank() == 0 )
1148  {
1149  /* divide pts to start */
1150 
1151  numPts = (int*)malloc( sizeof( int ) * mbpc->proc_config().proc_size() );
1152  ptsPerProc = npts / mbpc->proc_config().proc_size();
1153  ptsAssigned = 0;
1154 
1155  for( i = 0; i < mbpc->proc_config().proc_size() - 1; i++ )
1156  {
1157  numPts[i] = ptsPerProc;
1158  ptsAssigned += ptsPerProc;
1159  }
1160 
1161  numPts[mbpc->proc_config().proc_size() - 1] = npts - ptsAssigned;
1162 
1163  mySize = numPts[mbpc->proc_config().proc_rank()];
1164  sendPts = pts + ( 3 * numPts[0] );
1165  sendIds = ids + numPts[0];
1166  sum = 0; // possible no adjacency sent
1167  if( !part_geom )
1168  {
1169  sendEdges = length + numPts[0];
1170 
1171  for( j = 0; j < numPts[0]; j++ )
1172  sum += length[j];
1173 
1174  sendNborId = adjs + sum;
1175 
1176  for( j = numPts[0]; j < npts; j++ )
1177  sum += length[j];
1178 
1179  nborProcs = (int*)malloc( sizeof( int ) * sum );
1180  }
1181  for( j = 0; j < sum; j++ )
1182  if( ( i = adjs[j] / ptsPerProc ) < mbpc->proc_config().proc_size() )
1183  nborProcs[j] = i;
1184  else
1185  nborProcs[j] = mbpc->proc_config().proc_size() - 1;
1186 
1187  sendProcs = nborProcs + ( sendNborId - adjs );
1188 
1189  for( i = 1; i < mbpc->proc_config().proc_size(); i++ )
1190  {
1191  MPI_Send( &numPts[i], 1, MPI_INT, i, 0x00, mbpc->comm() );
1192  MPI_Send( sendPts, 3 * numPts[i], MPI_DOUBLE, i, 0x01, mbpc->comm() );
1193  MPI_Send( sendIds, numPts[i], MPI_INT, i, 0x03, mbpc->comm() );
1194  MPI_Send( sendEdges, numPts[i], MPI_INT, i, 0x06, mbpc->comm() );
1195  sum = 0;
1196 
1197  for( j = 0; j < numPts[i]; j++ )
1198  sum += sendEdges[j];
1199 
1200  MPI_Send( sendNborId, sum, MPI_INT, i, 0x07, mbpc->comm() );
1201  MPI_Send( sendProcs, sum, MPI_INT, i, 0x08, mbpc->comm() );
1202  sendPts += ( 3 * numPts[i] );
1203  sendIds += numPts[i];
1204  sendEdges += numPts[i];
1205  sendNborId += sum;
1206  sendProcs += sum;
1207  }
1208 
1209  free( numPts );
1210  }
1211  else
1212  {
1213  MPI_Recv( &mySize, 1, MPI_INT, 0, 0x00, mbpc->comm(), &stat );
1214  pts = (double*)malloc( sizeof( double ) * 3 * mySize );
1215  ids = (int*)malloc( sizeof( int ) * mySize );
1216  length = (int*)malloc( sizeof( int ) * mySize );
1217  if( obj_weights != NULL ) obj_weights = (double*)malloc( sizeof( double ) * mySize );
1218  MPI_Recv( pts, 3 * mySize, MPI_DOUBLE, 0, 0x01, mbpc->comm(), &stat );
1219  MPI_Recv( ids, mySize, MPI_INT, 0, 0x03, mbpc->comm(), &stat );
1220  MPI_Recv( length, mySize, MPI_INT, 0, 0x06, mbpc->comm(), &stat );
1221  sum = 0;
1222 
1223  for( j = 0; j < mySize; j++ )
1224  sum += length[j];
1225 
1226  adjs = (int*)malloc( sizeof( int ) * sum );
1227  if( edge_weights != NULL ) edge_weights = (double*)malloc( sizeof( double ) * sum );
1228  nborProcs = (int*)malloc( sizeof( int ) * sum );
1229  MPI_Recv( adjs, sum, MPI_INT, 0, 0x07, mbpc->comm(), &stat );
1230  MPI_Recv( nborProcs, sum, MPI_INT, 0, 0x08, mbpc->comm(), &stat );
1231  }
1232 
1233  Points = pts;
1234  GlobalIds = ids;
1235  NumPoints = mySize;
1236  NumEdges = length;
1237  NborGlobalId = adjs;
1238  NborProcs = nborProcs;
1239  ObjWeights = obj_weights;
1240  EdgeWeights = edge_weights;
1241  Parts = parts;
1242 
1243  return mySize;
1244 }
1245 
1246 void ZoltanPartitioner::mbFinalizePoints( int npts,
1247  int numExport,
1248  ZOLTAN_ID_PTR exportLocalIDs,
1249  int* exportProcs,
1250  int** assignment )
1251 {
1252  int* MyAssignment;
1253  int i;
1254  int numPts;
1255  MPI_Status stat;
1256  int* recvA;
1257 
1258  /* assign pts to start */
1259 
1260  if( mbpc->proc_config().proc_rank() == 0 )
1261  MyAssignment = (int*)malloc( sizeof( int ) * npts );
1262  else
1263  MyAssignment = (int*)malloc( sizeof( int ) * NumPoints );
1264 
1265  for( i = 0; i < NumPoints; i++ )
1266  MyAssignment[i] = mbpc->proc_config().proc_rank();
1267 
1268  for( i = 0; i < numExport; i++ )
1269  MyAssignment[exportLocalIDs[i]] = exportProcs[i];
1270 
1271  if( mbpc->proc_config().proc_rank() == 0 )
1272  {
1273  /* collect pts */
1274  recvA = MyAssignment + NumPoints;
1275 
1276  for( i = 1; i < (int)mbpc->proc_config().proc_size(); i++ )
1277  {
1278  MPI_Recv( &numPts, 1, MPI_INT, i, 0x04, mbpc->comm(), &stat );
1279  MPI_Recv( recvA, numPts, MPI_INT, i, 0x05, mbpc->comm(), &stat );
1280  recvA += numPts;
1281  }
1282 
1283  *assignment = MyAssignment;
1284  }
1285  else
1286  {
1287  MPI_Send( &NumPoints, 1, MPI_INT, 0, 0x04, mbpc->comm() );
1288  MPI_Send( MyAssignment, NumPoints, MPI_INT, 0, 0x05, mbpc->comm() );
1289  free( MyAssignment );
1290  }
1291 }
1292 
1293 int ZoltanPartitioner::mbGlobalSuccess( int rc )
1294 {
1295  int fail = 0;
1296  unsigned int i;
1297  int* vals = (int*)malloc( mbpc->proc_config().proc_size() * sizeof( int ) );
1298 
1299  MPI_Allgather( &rc, 1, MPI_INT, vals, 1, MPI_INT, mbpc->comm() );
1300 
1301  for( i = 0; i < mbpc->proc_config().proc_size(); i++ )
1302  {
1303  if( vals[i] != ZOLTAN_OK )
1304  {
1305  if( 0 == mbpc->proc_config().proc_rank() )
1306  {
1307  mbShowError( vals[i], "Result on process " );
1308  }
1309  fail = 1;
1310  }
1311  }
1312 
1313  free( vals );
1314  return fail;
1315 }
1316 
1317 void ZoltanPartitioner::mbPrintGlobalResult( const char* s, int begin, int import, int exp, int change )
1318 {
1319  unsigned int i;
1320  int* v1 = (int*)malloc( 4 * sizeof( int ) );
1321  int* v2 = NULL;
1322  int* v;
1323 
1324  v1[0] = begin;
1325  v1[1] = import;
1326  v1[2] = exp;
1327  v1[3] = change;
1328 
1329  if( mbpc->proc_config().proc_rank() == 0 )
1330  {
1331  v2 = (int*)malloc( 4 * mbpc->proc_config().proc_size() * sizeof( int ) );
1332  }
1333 
1334  MPI_Gather( v1, 4, MPI_INT, v2, 4, MPI_INT, 0, mbpc->comm() );
1335 
1336  if( mbpc->proc_config().proc_rank() == 0 )
1337  {
1338  fprintf( stdout, "======%s======\n", s );
1339  for( i = 0, v = v2; i < mbpc->proc_config().proc_size(); i++, v += 4 )
1340  {
1341  fprintf( stdout, "%u: originally had %d, import %d, exp %d, %s\n", i, v[0], v[1], v[2],
1342  v[3] ? "a change of partitioning" : "no change" );
1343  }
1344  fprintf( stdout, "==========================================\n" );
1345  fflush( stdout );
1346 
1347  free( v2 );
1348  }
1349 
1350  free( v1 );
1351 }
1352 
1353 void ZoltanPartitioner::mbShowError( int val, const char* s )
1354 {
1355  if( s ) printf( "%s ", s );
1356 
1357  switch( val )
1358  {
1359  case ZOLTAN_OK:
1360  printf( "%d: SUCCESSFUL\n", mbpc->proc_config().proc_rank() );
1361  break;
1362  case ZOLTAN_WARN:
1363  printf( "%d: WARNING\n", mbpc->proc_config().proc_rank() );
1364  break;
1365  case ZOLTAN_FATAL:
1366  printf( "%d: FATAL ERROR\n", mbpc->proc_config().proc_rank() );
1367  break;
1368  case ZOLTAN_MEMERR:
1369  printf( "%d: MEMORY ALLOCATION ERROR\n", mbpc->proc_config().proc_rank() );
1370  break;
1371  default:
1372  printf( "%d: INVALID RETURN CODE\n", mbpc->proc_config().proc_rank() );
1373  break;
1374  }
1375 }
1376 
1377 /**********************
1378 ** call backs
1379 **********************/
1380 
1381 int mbGetNumberOfAssignedObjects( void* /* userDefinedData */, int* err )
1382 {
1383  *err = 0;
1384  return NumPoints;
1385 }
1386 
1387 void mbGetObjectList( void* /* userDefinedData */,
1388  int /* numGlobalIds */,
1389  int /* numLids */,
1390  ZOLTAN_ID_PTR gids,
1391  ZOLTAN_ID_PTR lids,
1392  int wgt_dim,
1393  float* obj_wgts,
1394  int* err )
1395 {
1396  for( int i = 0; i < NumPoints; i++ )
1397  {
1398  gids[i] = GlobalIds[i];
1399  lids[i] = i;
1400  if( wgt_dim > 0 ) obj_wgts[i] = ObjWeights[i];
1401  }
1402 
1403  *err = 0;
1404 }
1405 
1406 int mbGetObjectSize( void* /* userDefinedData */, int* err )
1407 {
1408  *err = 0;
1409  return 3;
1410 }
1411 
1412 void mbGetObject( void* /* userDefinedData */,
1413  int /* numGlobalIds */,
1414  int /* numLids */,
1415  int numObjs,
1416  ZOLTAN_ID_PTR /* gids */,
1417  ZOLTAN_ID_PTR lids,
1418  int numDim,
1419  double* pts,
1420  int* err )
1421 {
1422  int i, id, id3;
1423  int next = 0;
1424 
1425  if( numDim != 3 )
1426  {
1427  *err = 1;
1428  return;
1429  }
1430 
1431  for( i = 0; i < numObjs; i++ )
1432  {
1433  id = lids[i];
1434 
1435  if( ( id < 0 ) || ( id >= NumPoints ) )
1436  {
1437  *err = 1;
1438  return;
1439  }
1440 
1441  id3 = lids[i] * 3;
1442 
1443  pts[next++] = Points[id3];
1444  pts[next++] = Points[id3 + 1];
1445  pts[next++] = Points[id3 + 2];
1446  }
1447 }
1448 
1449 void mbGetNumberOfEdges( void* /* userDefinedData */,
1450  int /* numGlobalIds */,
1451  int /* numLids */,
1452  int numObjs,
1453  ZOLTAN_ID_PTR /* gids */,
1454  ZOLTAN_ID_PTR lids,
1455  int* numEdges,
1456  int* err )
1457 {
1458  int i, id;
1459  int next = 0;
1460 
1461  for( i = 0; i < numObjs; i++ )
1462  {
1463  id = lids[i];
1464 
1465  if( ( id < 0 ) || ( id >= NumPoints ) )
1466  {
1467  *err = 1;
1468  return;
1469  }
1470 
1471  numEdges[next++] = NumEdges[id];
1472  }
1473 }
1474 
1475 void mbGetEdgeList( void* /* userDefinedData */,
1476  int /* numGlobalIds */,
1477  int /* numLids */,
1478  int numObjs,
1479  ZOLTAN_ID_PTR /* gids */,
1480  ZOLTAN_ID_PTR lids,
1481  int* /* numEdges */,
1482  ZOLTAN_ID_PTR nborGlobalIds,
1483  int* nborProcs,
1484  int wgt_dim,
1485  float* edge_wgts,
1486  int* err )
1487 {
1488  int i, id, idSum, j;
1489  int next = 0;
1490 
1491  for( i = 0; i < numObjs; i++ )
1492  {
1493  id = lids[i];
1494 
1495  if( ( id < 0 ) || ( id >= NumPoints ) )
1496  {
1497  *err = 1;
1498  return;
1499  }
1500 
1501  idSum = 0;
1502 
1503  for( j = 0; j < id; j++ )
1504  idSum += NumEdges[j];
1505 
1506  for( j = 0; j < NumEdges[id]; j++ )
1507  {
1508  nborGlobalIds[next] = NborGlobalId[idSum];
1509  nborProcs[next] = NborProcs[idSum];
1510  if( wgt_dim > 0 ) edge_wgts[next] = EdgeWeights[idSum];
1511  next++;
1512  idSum++;
1513  }
1514  }
1515 }
1516 
1517 void mbGetPart( void* /* userDefinedData */,
1518  int /* numGlobalIds */,
1519  int /* numLids */,
1520  int numObjs,
1521  ZOLTAN_ID_PTR /* gids */,
1522  ZOLTAN_ID_PTR lids,
1523  int* part,
1524  int* err )
1525 {
1526  int i, id;
1527  int next = 0;
1528 
1529  for( i = 0; i < numObjs; i++ )
1530  {
1531  id = lids[i];
1532 
1533  if( ( id < 0 ) || ( id >= NumPoints ) )
1534  {
1535  *err = 1;
1536  return;
1537  }
1538 
1539  part[next++] = Parts[id];
1540  }
1541 }
1542 
1543 // new methods for partition in parallel, used by migrate in iMOAB
1544 ErrorCode ZoltanPartitioner::partition_owned_cells( Range& primary,
1545  std::multimap< int, int >& extraGraphEdges,
1546  std::map< int, int > procs,
1547  int& numNewPartitions,
1548  std::map< int, Range >& distribution,
1549  int partition_method )
1550 {
1551  // start copy
1552  MeshTopoUtil mtu( mbImpl );
1553  Range adjs;
1554 
1555  std::vector< int > adjacencies;
1556  std::vector< int > ids;
1557  ids.resize( primary.size() );
1558  std::vector< int > length;
1559  std::vector< double > coords;
1560  std::vector< int > nbor_proc;
1561  // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited
1562  // by MBCN
1563  int neighbors[5 * MAX_SUB_ENTITIES];
1564  int neib_proc[5 * MAX_SUB_ENTITIES];
1565  double avg_position[3];
1566  int moab_id;
1567  int primaryDim = mbImpl->dimension_from_handle( *primary.rbegin() );
1568  // get the global id tag handle
1569  Tag gid = mbImpl->globalId_tag();
1570 
1571  MB_CHK_ERR( mbImpl->tag_get_data( gid, primary, &ids[0] ) );
1572 
1573  // mbpc is member in base class, PartitionerBase
1574  int rank = mbpc->rank(); // current rank , will be put on regular neighbors
1575  int i = 0;
1576  for( Range::iterator rit = primary.begin(); rit != primary.end(); ++rit, i++ )
1577  {
1578  EntityHandle cell = *rit;
1579  // get bridge adjacencies for each cell
1580  if( 1 == partition_method )
1581  {
1582  adjs.clear();
1583  MB_CHK_ERR( mtu.get_bridge_adjacencies( cell, ( primaryDim > 0 ? primaryDim - 1 : 3 ), primaryDim, adjs ) );
1584 
1585  // get the graph vertex ids of those
1586  if( !adjs.empty() )
1587  {
1588  assert( adjs.size() < 5 * MAX_SUB_ENTITIES );
1589  MB_CHK_ERR( mbImpl->tag_get_data( gid, adjs, neighbors ) );
1590  }
1591  // if adjacent to neighbor partitions, add to the list
1592  int size_adjs = (int)adjs.size();
1593  moab_id = ids[i];
1594  for( int k = 0; k < size_adjs; k++ )
1595  neib_proc[k] = rank; // current rank
1596  if( extraGraphEdges.find( moab_id ) != extraGraphEdges.end() )
1597  {
1598  // it means that the current cell is adjacent to a cell in another partition ; maybe
1599  // a few
1600  std::pair< std::multimap< int, int >::iterator, std::multimap< int, int >::iterator > ret;
1601  ret = extraGraphEdges.equal_range( moab_id );
1602  for( std::multimap< int, int >::iterator it = ret.first; it != ret.second; ++it )
1603  {
1604  int otherID = it->second;
1605  neighbors[size_adjs] = otherID; // the id of the other cell, across partition
1606  neib_proc[size_adjs] = procs[otherID]; // this is how we built this map, the cell id maps to
1607  // what proc it came from
1608  size_adjs++;
1609  }
1610  }
1611  // copy those into adjacencies vector
1612  length.push_back( size_adjs );
1613  std::copy( neighbors, neighbors + size_adjs, std::back_inserter( adjacencies ) );
1614  std::copy( neib_proc, neib_proc + size_adjs, std::back_inserter( nbor_proc ) );
1615  }
1616  else if( 2 == partition_method )
1617  {
1618  if( TYPE_FROM_HANDLE( cell ) == MBVERTEX )
1619  {
1620  MB_CHK_ERR( mbImpl->get_coords( &cell, 1, avg_position ) );
1621  }
1622  else
1623  {
1624  MB_CHK_ERR( mtu.get_average_position( cell, avg_position ) );
1625  }
1626  std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) );
1627  }
1628  }
1629 
1630 #ifdef VERBOSE
1631  std::stringstream ff2;
1632  ff2 << "zoltanInput_" << mbpc->rank() << ".txt";
1633  std::ofstream ofs;
1634  ofs.open( ff2.str().c_str(), std::ofstream::out );
1635  ofs << "Length vector: " << std::endl;
1636  std::copy( length.begin(), length.end(), std::ostream_iterator< int >( ofs, ", " ) );
1637  ofs << std::endl;
1638  ofs << "Adjacencies vector: " << std::endl;
1639  std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< int >( ofs, ", " ) );
1640  ofs << std::endl;
1641  ofs << "Moab_ids vector: " << std::endl;
1642  std::copy( ids.begin(), ids.end(), std::ostream_iterator< int >( ofs, ", " ) );
1643  ofs << std::endl;
1644  ofs << "Coords vector: " << std::endl;
1645  std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( ofs, ", " ) );
1646  ofs << std::endl;
1647  ofs.close();
1648 #endif
1649 
1650  // these are static var in this file, and used in the callbacks
1651  Points = NULL;
1652  if( 1 != partition_method ) Points = &coords[0];
1653  GlobalIds = &ids[0];
1654  NumPoints = (int)ids.size();
1655  NumEdges = &length[0];
1656  NborGlobalId = &adjacencies[0];
1657  NborProcs = &nbor_proc[0];
1658  ObjWeights = NULL;
1659  EdgeWeights = NULL;
1660  Parts = NULL;
1661 
1662  float version;
1663  Zoltan_Initialize( argcArg, argvArg, &version );
1664 
1665  // Create Zoltan object. This calls Zoltan_Create.
1666  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
1667 
1668  // set # requested partitions
1669  char buff[10];
1670  snprintf( buff, 10, "%d", numNewPartitions );
1671  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
1672  if( ZOLTAN_OK != retval ) return MB_FAILURE;
1673 
1674  // request parts assignment
1675  retval = myZZ->Set_Param( "RETURN_LISTS", "PARTS" );
1676  if( ZOLTAN_OK != retval ) return MB_FAILURE;
1677 
1678  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
1679  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
1680  // due to a bug in zoltan, if method is graph partitioning, do not pass coordinates!!
1681  if( 2 == partition_method )
1682  {
1683  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
1684  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
1685  SetRCB_Parameters(); // geometry
1686  }
1687  else if( 1 == partition_method )
1688  {
1689  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
1690  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
1691  SetHypergraph_Parameters( "auto" );
1692  }
1693 
1694  // Perform the load balancing partitioning
1695 
1696  int changes;
1697  int numGidEntries;
1698  int numLidEntries;
1699  int num_import;
1700  ZOLTAN_ID_PTR import_global_ids, import_local_ids;
1701  int* import_procs;
1702  int* import_to_part;
1703  int num_export;
1704  ZOLTAN_ID_PTR export_global_ids, export_local_ids;
1705  int *assign_procs, *assign_parts;
1706 
1707  // if( mbpc->rank() == 0 )
1708  // std::cout << "Computing partition using method (1-graph, 2-geom): " << partition_method << " for " << numNewPartitions
1709  // << " parts..." << std::endl;
1710 
1711 #ifndef NDEBUG
1712 #if 0
1713  static int counter=0; // it may be possible to call function multiple times in a simulation
1714  // give a way to not overwrite the files
1715  // it should work only with a modified version of Zoltan
1716  std::stringstream basename;
1717  if (1==partition_method)
1718  {
1719  basename << "phg_" << counter++;
1720  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)(basename.str().c_str()), 1, 0, 1, 0);
1721  }
1722  else if (2==partition_method)
1723  {
1724  basename << "rcb_" << counter++;
1725  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)(basename.str().c_str()), 1, 1, 0, 0);
1726  }
1727 #endif
1728 #endif
1729  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids,
1730  import_procs, import_to_part, num_export, export_global_ids, export_local_ids,
1731  assign_procs, assign_parts );
1732  if( ZOLTAN_OK != retval ) return MB_FAILURE;
1733 
1734 #ifdef VERBOSE
1735  std::stringstream ff3;
1736  ff3 << "zoltanOutput_" << mbpc->rank() << ".txt";
1737  std::ofstream ofs3;
1738  ofs3.open( ff3.str().c_str(), std::ofstream::out );
1739  ofs3 << " export elements on rank " << rank << " \n";
1740  ofs3 << "\t index \t gb_id \t local \t proc \t part \n";
1741  for( int k = 0; k < num_export; k++ )
1742  {
1743  ofs3 << "\t" << k << "\t" << export_global_ids[k] << "\t" << export_local_ids[k] << "\t" << assign_procs[k]
1744  << "\t" << assign_parts[k] << "\n";
1745  }
1746  ofs3.close();
1747 #endif
1748 
1749  // basically each local cell is assigned to a part
1750 
1751  assert( num_export == (int)primary.size() );
1752  for( i = 0; i < num_export; i++ )
1753  {
1754  EntityHandle cell = primary[export_local_ids[i]];
1755  distribution[assign_parts[i]].insert( cell );
1756  }
1757 
1758  Zoltan::LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part );
1759  Zoltan::LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts );
1760 
1761  delete myZZ;
1762  myZZ = NULL;
1763 
1764  // clear arrays that were resized locally, to free up local memory
1765 
1766  std::vector< int >().swap( adjacencies );
1767  std::vector< int >().swap( ids );
1768  std::vector< int >().swap( length );
1769  std::vector< int >().swap( nbor_proc );
1770  std::vector< double >().swap( coords );
1771 
1772  return MB_SUCCESS;
1773 }