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