EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_MapColoring.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 
42 #include <Epetra_ConfigDefs.h>
43 #include <EpetraExt_MapColoring.h>
44 
47 
48 #include <Epetra_CrsGraph.h>
49 #include <Epetra_GIDTypeVector.h>
50 #include <Epetra_MapColoring.h>
51 #include <Epetra_Map.h>
52 #include <Epetra_Comm.h>
53 #include <Epetra_Util.h>
54 #include <Epetra_Import.h>
55 #include <Epetra_Export.h>
56 
57 #include <Epetra_Time.h>
58 
59 #include <vector>
60 #include <set>
61 #include <map>
62 
63 using std::vector;
64 using std::set;
65 using std::map;
66 
67 namespace EpetraExt {
68 
69 template<typename int_type>
71 CrsGraph_MapColoring::
72 Toperator( OriginalTypeRef orig )
73 {
74  Epetra_Time timer( orig.Comm() );
75 
76  origObj_ = &orig;
77 
78  const Epetra_BlockMap & RowMap = orig.RowMap();
79  int nRows = RowMap.NumMyElements();
80  const Epetra_BlockMap & ColMap = orig.ColMap();
81  int nCols = ColMap.NumMyElements();
82 
83  int MyPID = RowMap.Comm().MyPID();
84 
85  if( verbosity_ > 1 ) std::cout << "RowMap:\n" << RowMap;
86  if( verbosity_ > 1 ) std::cout << "ColMap:\n" << ColMap;
87 
88  Epetra_MapColoring * ColorMap = 0;
89 
90  if( algo_ == GREEDY || algo_ == LUBY )
91  {
92 
93  Epetra_CrsGraph * base = &( const_cast<Epetra_CrsGraph&>(orig) );
94 
95  int NumIndices;
96  int * Indices;
97 
98  double wTime1 = timer.WallTime();
99 
100  // For parallel applications, add in boundaries to coloring
101  bool distributedGraph = RowMap.DistributedGlobal();
102  if( distributedGraph )
103  {
104  base = new Epetra_CrsGraph( Copy, ColMap, ColMap, 0 );
105 
106  for( int i = 0; i < nRows; ++i )
107  {
108  orig.ExtractMyRowView( i, NumIndices, Indices );
109  base->InsertMyIndices( i, NumIndices, Indices );
110 
111  //Do this with a single insert
112  //Is this the right thing?
113  for( int j = 0; j < NumIndices; ++j )
114  if( Indices[j] >= nRows )
115  base->InsertMyIndices( Indices[j], 1, &i );
116  }
117  base->FillComplete();
118  }
119 
120  if( verbosity_ > 1 ) std::cout << "Base Graph:\n" << *base << std::endl;
121 
122  double wTime2 = timer.WallTime();
123  if( verbosity_ > 0 )
124  std::cout << "EpetraExt::MapColoring [INSERT BOUNDARIES] Time: " << wTime2-wTime1 << std::endl;
125 
126  //Generate Local Distance-1 Adjacency Graph
127  //(Transpose of orig excluding non-local column indices)
128  EpetraExt::CrsGraph_Transpose transposeTransform( true );
129  Epetra_CrsGraph & Adj1 = transposeTransform( *base );
130  if( verbosity_ > 1 ) std::cout << "Adjacency 1 Graph!\n" << Adj1;
131 
132  wTime1 = timer.WallTime();
133  if( verbosity_ > 0 )
134  std::cout << "EpetraExt::MapColoring [TRANSPOSE GRAPH] Time: " << wTime1-wTime2 << std::endl;
135 
136  int Delta = Adj1.MaxNumIndices();
137  if( verbosity_ > 0 ) std::cout << std::endl << "Delta: " << Delta << std::endl;
138 
139  //Generation of Local Distance-2 Adjacency Graph
140  Epetra_CrsGraph * Adj2 = &Adj1;
141  if( !distance1_ )
142  {
143  Adj2 = new Epetra_CrsGraph( Copy, ColMap, ColMap, 0 );
144  int NumAdj1Indices;
145  int * Adj1Indices;
146  for( int i = 0; i < nCols; ++i )
147  {
148  Adj1.ExtractMyRowView( i, NumAdj1Indices, Adj1Indices );
149 
150  set<int> Cols;
151  for( int j = 0; j < NumAdj1Indices; ++j )
152  {
153  base->ExtractMyRowView( Adj1Indices[j], NumIndices, Indices );
154  for( int k = 0; k < NumIndices; ++k )
155  if( Indices[k] < nCols ) Cols.insert( Indices[k] );
156  }
157  int nCols2 = Cols.size();
158  std::vector<int> ColVec( nCols2 );
159  set<int>::iterator iterIS = Cols.begin();
160  set<int>::iterator iendIS = Cols.end();
161  for( int j = 0 ; iterIS != iendIS; ++iterIS, ++j ) ColVec[j] = *iterIS;
162  Adj2->InsertMyIndices( i, nCols2, &ColVec[0] );
163  }
164  Adj2->FillComplete();
165 
166  if( verbosity_ > 1 ) std::cout << "Adjacency 2 Graph!\n" << *Adj2;
167 
168  wTime2 = timer.WallTime();
169  if( verbosity_ > 0 )
170  std::cout << "EpetraExt::MapColoring [GEN DIST-2 GRAPH] Time: " << wTime2-wTime1 << std::endl;
171  }
172 
173  wTime2 = timer.WallTime();
174 
175  ColorMap = new Epetra_MapColoring( ColMap );
176 
177  std::vector<int> rowOrder( nCols );
178  if( reordering_ == 0 || reordering_ == 1 )
179  {
180  std::multimap<int,int> adjMap;
181  typedef std::multimap<int,int>::value_type adjMapValueType;
182  for( int i = 0; i < nCols; ++i )
183  adjMap.insert( adjMapValueType( Adj2->NumMyIndices(i), i ) );
184  std::multimap<int,int>::iterator iter = adjMap.begin();
185  std::multimap<int,int>::iterator end = adjMap.end();
186  if( reordering_ == 0 ) //largest first (less colors)
187  {
188  for( int i = 1; iter != end; ++iter, ++i )
189  rowOrder[ nCols - i ] = (*iter).second;
190  }
191  else //smallest first (better balance)
192  {
193  for( int i = 0; iter != end; ++iter, ++i )
194  rowOrder[ i ] = (*iter).second;
195  }
196  }
197  else if( reordering_ == 2 ) //random
198  {
199  for( int i = 0; i < nCols; ++i )
200  rowOrder[ i ] = i;
201 #ifndef TFLOP
202  std::random_shuffle( rowOrder.begin(), rowOrder.end() );
203 #endif
204  }
205 
206  wTime1 = timer.WallTime();
207  if( verbosity_ > 0 )
208  std::cout << "EpetraExt::MapColoring [REORDERING] Time: " << wTime1-wTime2 << std::endl;
209 
210  //Application of Greedy Algorithm to generate Color Map
211  if( algo_ == GREEDY )
212  {
213  for( int col = 0; col < nCols; ++col )
214  {
215  Adj2->ExtractMyRowView( rowOrder[col], NumIndices, Indices );
216 
217  set<int> usedColors;
218  int color;
219  for( int i = 0; i < NumIndices; ++i )
220  {
221  color = (*ColorMap)[ Indices[i] ];
222  if( color > 0 ) usedColors.insert( color );
223  color = 0;
224  int testcolor = 1;
225  while( !color )
226  {
227  if( !usedColors.count( testcolor ) ) color = testcolor;
228  ++testcolor;
229  }
230  }
231  (*ColorMap)[ rowOrder[col] ] = color;
232  }
233 
234  wTime2 = timer.WallTime();
235  if( verbosity_ > 0 )
236  std::cout << "EpetraExt::MapColoring [GREEDY COLORING] Time: " << wTime2-wTime1 << std::endl;
237  if( verbosity_ > 0 )
238  std::cout << "Num GREEDY Colors: " << ColorMap->NumColors() << std::endl;
239  }
240  else if( algo_ == LUBY )
241  {
242  //Assign Random Keys To Rows
243  Epetra_Util util;
244  std::vector<int> Keys(nCols);
245  std::vector<int> State(nCols,-1);
246 
247  for( int col = 0; col < nCols; ++col )
248  Keys[col] = util.RandomInt();
249 
250  int NumRemaining = nCols;
251  int CurrentColor = 1;
252 
253  while( NumRemaining > 0 )
254  {
255  //maximal independent set
256  while( NumRemaining > 0 )
257  {
258  NumRemaining = 0;
259 
260  //zero out everyone less than neighbor
261  for( int i = 0; i < nCols; ++i )
262  {
263  int col = rowOrder[i];
264  if( State[col] < 0 )
265  {
266  Adj2->ExtractMyRowView( col, NumIndices, Indices );
267  int MyKey = Keys[col];
268  for( int j = 0; j < NumIndices; ++j )
269  if( col != Indices[j] && State[ Indices[j] ] < 0 )
270  {
271  if( MyKey > Keys[ Indices[j] ] ) State[ Indices[j] ] = 0;
272  else State[ col ] = 0;
273  }
274  }
275  }
276 
277  //assign -1's to current color
278  for( int col = 0; col < nCols; ++col )
279  {
280  if( State[col] < 0 )
281  State[col] = CurrentColor;
282  }
283 
284  //reinstate any zero not neighboring current color
285  for( int col = 0; col < nCols; ++col )
286  if( State[col] == 0 )
287  {
288  Adj2->ExtractMyRowView( col, NumIndices, Indices );
289  bool flag = false;
290  for( int i = 0; i < NumIndices; ++i )
291  if( col != Indices[i] && State[ Indices[i] ] == CurrentColor )
292  {
293  flag = true;
294  break;
295  }
296  if( !flag )
297  {
298  State[col] = -1;
299  ++NumRemaining;
300  }
301  }
302  }
303 
304  //Reset Status for all non-colored nodes
305  for( int col = 0; col < nCols; ++col )
306  if( State[col] == 0 )
307  {
308  State[col] = -1;
309  ++NumRemaining;
310  }
311 
312  if( verbosity_ > 2 )
313  {
314  std::cout << "Finished Color: " << CurrentColor << std::endl;
315  std::cout << "NumRemaining: " << NumRemaining << std::endl;
316  }
317 
318  //New color
319  ++CurrentColor;
320  }
321 
322  for( int col = 0; col < nCols; ++col )
323  (*ColorMap)[col] = State[col]-1;
324 
325  wTime2 = timer.WallTime();
326  if( verbosity_ > 0 )
327  std::cout << "EpetraExt::MapColoring [LUBI COLORING] Time: " << wTime2-wTime1 << std::endl;
328  if( verbosity_ > 0 )
329  std::cout << "Num LUBI Colors: " << ColorMap->NumColors() << std::endl;
330 
331  }
332  else
333  abort(); //UNKNOWN ALGORITHM
334 
335  if( distributedGraph ) delete base;
336  if( !distance1_ ) delete Adj2;
337  }
338  else
339  {
340  //Generate Parallel Adjacency-2 Graph
341  EpetraExt::CrsGraph_Overlap OverlapTrans(1);
342  Epetra_CrsGraph & OverlapGraph = OverlapTrans( orig );
343 
344  if( verbosity_ > 1 ) std::cout << "OverlapGraph:\n" << OverlapGraph;
345 
346  Epetra_CrsGraph * Adj2 = &orig;
347 
348  int NumIndices;
349  int * Indices;
350  if( !distance1_ )
351  {
352  Adj2 = new Epetra_CrsGraph( Copy, RowMap, OverlapGraph.ColMap(), 0 );
353  int NumAdj1Indices;
354  int * Adj1Indices;
355  for( int i = 0; i < nRows; ++i )
356  {
357  OverlapGraph.ExtractMyRowView( i, NumAdj1Indices, Adj1Indices );
358 
359  set<int> Cols;
360  for( int j = 0; j < NumAdj1Indices; ++j )
361  {
362  int GID = OverlapGraph.LRID( OverlapGraph.GCID64( Adj1Indices[j] ) );
363  OverlapGraph.ExtractMyRowView( GID, NumIndices, Indices );
364  for( int k = 0; k < NumIndices; ++k ) Cols.insert( Indices[k] );
365  }
366  int nCols2 = Cols.size();
367  std::vector<int> ColVec( nCols2 );
368  set<int>::iterator iterIS = Cols.begin();
369  set<int>::iterator iendIS = Cols.end();
370  for( int j = 0 ; iterIS != iendIS; ++iterIS, ++j ) ColVec[j] = *iterIS;
371  Adj2->InsertMyIndices( i, nCols2, &ColVec[0] );
372  }
373 
374 #ifdef NDEBUG
375  (void) Adj2->FillComplete();
376 #else
377  // assert() statements go away if NDEBUG is defined. Don't
378  // declare the 'flag' variable if it never gets used.
379  int flag = Adj2->FillComplete();
380  assert( flag == 0 );
381 #endif // NDEBUG
382 
383  RowMap.Comm().Barrier();
384  if( verbosity_ > 1 ) std::cout << "Adjacency 2 Graph!\n" << *Adj2;
385  }
386 
387  //collect GIDs on boundary
388  std::vector<int_type> boundaryGIDs;
389  std::vector<int_type> interiorGIDs;
390  for( int row = 0; row < nRows; ++row )
391  {
392  Adj2->ExtractMyRowView( row, NumIndices, Indices );
393  bool testFlag = false;
394  for( int i = 0; i < NumIndices; ++i )
395  if( Indices[i] >= nRows ) testFlag = true;
396  if( testFlag ) boundaryGIDs.push_back( (int_type) Adj2->GRID64(row) );
397  else interiorGIDs.push_back( (int_type) Adj2->GRID64(row) );
398  }
399 
400  int LocalBoundarySize = boundaryGIDs.size();
401 
402  Epetra_Map BoundaryMap( -1, boundaryGIDs.size(),
403  LocalBoundarySize ? &boundaryGIDs[0]: 0,
404  0, RowMap.Comm() );
405  if( verbosity_ > 1 ) std::cout << "BoundaryMap:\n" << BoundaryMap;
406 
407  int_type BoundarySize = (int_type) BoundaryMap.NumGlobalElements64();
408  Epetra_MapColoring BoundaryColoring( BoundaryMap );
409 
410  if( algo_ == PSEUDO_PARALLEL )
411  {
412  Epetra_Map BoundaryIndexMap( BoundarySize, LocalBoundarySize, 0, RowMap.Comm() );
413  if( verbosity_ > 1) std::cout << "BoundaryIndexMap:\n" << BoundaryIndexMap;
414 
415  typename Epetra_GIDTypeVector<int_type>::impl bGIDs( View, BoundaryIndexMap, &boundaryGIDs[0] );
416  if( verbosity_ > 1) std::cout << "BoundaryGIDs:\n" << bGIDs;
417 
418  int_type NumLocalBs = 0;
419  if( !RowMap.Comm().MyPID() ) NumLocalBs = BoundarySize;
420 
421  Epetra_Map LocalBoundaryIndexMap( BoundarySize, NumLocalBs, 0, RowMap.Comm() );
422  if( verbosity_ > 1) std::cout << "LocalBoundaryIndexMap:\n" << LocalBoundaryIndexMap;
423 
424  typename Epetra_GIDTypeVector<int_type>::impl lbGIDs( LocalBoundaryIndexMap );
425  Epetra_Import lbImport( LocalBoundaryIndexMap, BoundaryIndexMap );
426  lbGIDs.Import( bGIDs, lbImport, Insert );
427  if( verbosity_ > 1) std::cout << "LocalBoundaryGIDs:\n" << lbGIDs;
428 
429  Epetra_Map LocalBoundaryMap( BoundarySize, NumLocalBs, lbGIDs.Values(), 0, RowMap.Comm() );
430  if( verbosity_ > 1) std::cout << "LocalBoundaryMap:\n" << LocalBoundaryMap;
431 
432  Epetra_CrsGraph LocalBoundaryGraph( Copy, LocalBoundaryMap, LocalBoundaryMap, 0 );
433  Epetra_Import LocalBoundaryImport( LocalBoundaryMap, Adj2->RowMap() );
434  LocalBoundaryGraph.Import( *Adj2, LocalBoundaryImport, Insert );
435  LocalBoundaryGraph.FillComplete();
436  if( verbosity_ > 1 ) std::cout << "LocalBoundaryGraph:\n " << LocalBoundaryGraph;
437 
438  EpetraExt::CrsGraph_MapColoring BoundaryTrans( GREEDY, reordering_, distance1_, verbosity_ );
439  Epetra_MapColoring & LocalBoundaryColoring = BoundaryTrans( LocalBoundaryGraph );
440  if( verbosity_ > 1 ) std::cout << "LocalBoundaryColoring:\n " << LocalBoundaryColoring;
441 
442  Epetra_Export BoundaryExport( LocalBoundaryMap, BoundaryMap );
443  BoundaryColoring.Export( LocalBoundaryColoring, BoundaryExport, Insert );
444  }
445  else if( algo_ == JONES_PLASSMAN )
446  {
447  /* Alternative Distrib. Memory Coloring of Boundary based on JonesPlassman(sic) paper
448  * 1.Random number assignment to all boundary nodes using GID as seed to function
449  * (This allows any processor to compute adj. off proc values with a local computation)
450  * 2.Find all nodes greater than any neighbor off processor, color them.
451  * 3.Send colored node info to neighbors
452  * 4.Constrained color all nodes with all off proc neighbors smaller or colored.
453  * 5.Goto 3
454  */
455 
456  std::vector<int_type> OverlapBoundaryGIDs( boundaryGIDs );
457  for( int i = nRows; i < Adj2->ColMap().NumMyElements(); ++i )
458  OverlapBoundaryGIDs.push_back( (int_type) Adj2->ColMap().GID64(i) );
459 
460  int_type OverlapBoundarySize = OverlapBoundaryGIDs.size();
461  Epetra_Map BoundaryColMap( (int_type) -1, OverlapBoundarySize,
462  OverlapBoundarySize ? &OverlapBoundaryGIDs[0] : 0,
463  0, RowMap.Comm() );
464 
465  Epetra_CrsGraph BoundaryGraph( Copy, BoundaryMap, BoundaryColMap, 0 );
466  Epetra_Import BoundaryImport( BoundaryMap, Adj2->RowMap() );
467  BoundaryGraph.Import( *Adj2, BoundaryImport, Insert );
468  BoundaryGraph.FillComplete();
469  if( verbosity_ > 1) std::cout << "BoundaryGraph:\n" << BoundaryGraph;
470 
471  Epetra_Import ReverseOverlapBoundaryImport( BoundaryMap, BoundaryColMap );
472  Epetra_Import OverlapBoundaryImport( BoundaryColMap, BoundaryMap );
473 
474  int Colored = 0;
475  int GlobalColored = 0;
476  int Level = 0;
477  Epetra_MapColoring OverlapBoundaryColoring( BoundaryColMap );
478 
479  //Setup random integers for boundary nodes
480  Epetra_IntVector BoundaryValues( BoundaryMap );
481  Epetra_Util Util;
482  Util.SetSeed( 47954118 * (MyPID+1) );
483  for( int i=0; i < LocalBoundarySize; ++i )
484  {
485  int val = Util.RandomInt();
486  if( val < 0 ) val *= -1;
487  BoundaryValues[i] = val;
488  }
489 
490  //Get Random Values for External Boundary
491  Epetra_IntVector OverlapBoundaryValues( BoundaryColMap );
492  OverlapBoundaryValues.Import( BoundaryValues, OverlapBoundaryImport, Insert );
493 
494  while( GlobalColored < BoundarySize )
495  {
496  //Find current "Level" of boundary indices to color
497  int theNumIndices;
498  int * theIndices;
499  std::vector<int> LevelIndices;
500  for( int i = 0; i < LocalBoundarySize; ++i )
501  {
502  if( !OverlapBoundaryColoring[i] )
503  {
504  //int MyVal = PRAND(BoundaryColMap.GID(i));
505  int MyVal = OverlapBoundaryValues[i];
506  BoundaryGraph.ExtractMyRowView( i, theNumIndices, theIndices );
507  bool ColorFlag = true;
508  int Loc = 0;
509  while( Loc<theNumIndices && theIndices[Loc]<LocalBoundarySize ) ++Loc;
510  for( int j = Loc; j < theNumIndices; ++j )
511  if( (OverlapBoundaryValues[theIndices[j]]>MyVal)
512  && !OverlapBoundaryColoring[theIndices[j]] )
513  {
514  ColorFlag = false;
515  break;
516  }
517  if( ColorFlag ) LevelIndices.push_back(i);
518  }
519  }
520 
521  if( verbosity_ > 1 )
522  {
523  std::cout << MyPID << " Level Indices: ";
524  int Lsize = (int) LevelIndices.size();
525  for( int i = 0; i < Lsize; ++i ) std::cout << LevelIndices[i] << " ";
526  std::cout << std::endl;
527  }
528 
529  //Greedy coloring of current level
530  set<int> levelColors;
531  int Lsize = (int) LevelIndices.size();
532  for( int i = 0; i < Lsize; ++i )
533  {
534  BoundaryGraph.ExtractMyRowView( LevelIndices[i], theNumIndices, theIndices );
535 
536  set<int> usedColors;
537  int color;
538  for( int j = 0; j < theNumIndices; ++j )
539  {
540  color = OverlapBoundaryColoring[ theIndices[j] ];
541  if( color > 0 ) usedColors.insert( color );
542  color = 0;
543  int testcolor = 1;
544  while( !color )
545  {
546  if( !usedColors.count( testcolor ) ) color = testcolor;
547  ++testcolor;
548  }
549  }
550  OverlapBoundaryColoring[ LevelIndices[i] ] = color;
551  levelColors.insert( color );
552  }
553 
554  if( verbosity_ > 2 ) std::cout << MyPID << " Level: " << Level << " Count: " << LevelIndices.size() << " NumColors: " << levelColors.size() << std::endl;
555 
556  if( verbosity_ > 2 ) std::cout << "Current Level Boundary Coloring:\n" << OverlapBoundaryColoring;
557 
558  //Update off processor coloring info
559  BoundaryColoring.Import( OverlapBoundaryColoring, ReverseOverlapBoundaryImport, Insert );
560  OverlapBoundaryColoring.Import( BoundaryColoring, OverlapBoundaryImport, Insert );
561  Colored += LevelIndices.size();
562  Level++;
563 
564  RowMap.Comm().SumAll( &Colored, &GlobalColored, 1 );
565  if( verbosity_ > 2 ) std::cout << "Num Globally Colored: " << GlobalColored << " from Num Global Boundary Nodes: " << BoundarySize << std::endl;
566  }
567  }
568 
569  if( verbosity_ > 1 ) std::cout << "BoundaryColoring:\n " << BoundaryColoring;
570 
571  Epetra_MapColoring RowColorMap( RowMap );
572 
573  //Add Boundary Colors
574  for( int i = 0; i < LocalBoundarySize; ++i )
575  {
576  int_type GID = (int_type) BoundaryMap.GID64(i);
577  RowColorMap(GID) = BoundaryColoring(GID);
578  }
579 
580  Epetra_MapColoring Adj2ColColorMap( Adj2->ColMap() );
581  Epetra_Import Adj2Import( Adj2->ColMap(), RowMap );
582  Adj2ColColorMap.Import( RowColorMap, Adj2Import, Insert );
583 
584  if( verbosity_ > 1 ) std::cout << "RowColoringMap:\n " << RowColorMap;
585  if( verbosity_ > 1 ) std::cout << "Adj2ColColorMap:\n " << Adj2ColColorMap;
586 
587  std::vector<int> rowOrder( nRows );
588  if( reordering_ == 0 || reordering_ == 1 )
589  {
590  std::multimap<int,int> adjMap;
591  typedef std::multimap<int,int>::value_type adjMapValueType;
592  for( int i = 0; i < nRows; ++i )
593  adjMap.insert( adjMapValueType( Adj2->NumMyIndices(i), i ) );
594  std::multimap<int,int>::iterator iter = adjMap.begin();
595  std::multimap<int,int>::iterator end = adjMap.end();
596  if( reordering_ == 0 ) //largest first (less colors)
597  {
598  for( int i = 1; iter != end; ++iter, ++i )
599  rowOrder[nRows-i] = (*iter).second;
600  }
601  else //smallest first (better balance)
602  {
603  for( int i = 0; iter != end; ++iter, ++i )
604  rowOrder[i] = (*iter).second;
605  }
606  }
607  else if( reordering_ == 2 ) //random
608  {
609  for( int i = 0; i < nRows; ++i )
610  rowOrder[i] = i;
611 #ifdef TFLOP
612  random_shuffle( rowOrder.begin(), rowOrder.end() );
613 #endif
614  }
615 
616  //Constrained greedy coloring of interior
617  set<int> InteriorColors;
618  for( int row = 0; row < nRows; ++row )
619  {
620  if( !RowColorMap[ rowOrder[row] ] )
621  {
622  Adj2->ExtractMyRowView( rowOrder[row], NumIndices, Indices );
623 
624  set<int> usedColors;
625  int color;
626  for( int i = 0; i < NumIndices; ++i )
627  {
628  color = Adj2ColColorMap[ Indices[i] ];
629  if( color > 0 ) usedColors.insert( color );
630  color = 0;
631  int testcolor = 1;
632  while( !color )
633  {
634  if( !usedColors.count( testcolor ) ) color = testcolor;
635  ++testcolor;
636  }
637  }
638  Adj2ColColorMap[ rowOrder[row] ] = color;
639  InteriorColors.insert( color );
640  }
641  }
642  if( verbosity_ > 1 ) std::cout << MyPID << " Num Interior Colors: " << InteriorColors.size() << std::endl;
643  if( verbosity_ > 1 ) std::cout << "RowColorMap after Greedy:\n " << RowColorMap;
644 
645  ColorMap = new Epetra_MapColoring( ColMap );
646  Epetra_Import ColImport( ColMap, Adj2->ColMap() );
647  ColorMap->Import( Adj2ColColorMap, ColImport, Insert );
648 
649  if( !distance1_ ) delete Adj2;
650  }
651 
652  if( verbosity_ > 0 ) std::cout << MyPID << " ColorMap Color Count: " << ColorMap->NumColors() << std::endl;
653  if( verbosity_ > 1 ) std::cout << "ColorMap!\n" << *ColorMap;
654 
655  newObj_ = ColorMap;
656 
657  return *ColorMap;
658 }
659 
663 {
664 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
665  if(orig.RowMap().GlobalIndicesInt()) {
666  return Toperator<int>(orig);
667  }
668  else
669 #endif
670 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
671  if(orig.RowMap().GlobalIndicesLongLong()) {
672  return Toperator<long long>(orig);
673  }
674  else
675 #endif
676  throw "CrsGraph_MapColoring::operator(): ERROR, GlobalIndices type unknown.";
677 }
678 
679 } // namespace EpetraExt
bool DistributedGlobal() const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
Transform to generate the explicit transpose of a Epetra_CrsGraph.
const Epetra_BlockMap & ColMap() const
Given an input Epetra_CrsGraph, a &quot;overlapped&quot; Epetra_CrsGraph is generated including rows associated...
virtual int MyPID() const =0
unsigned int RandomInt()
CrsGraph_MapColoring::NewTypeRef operator()(CrsGraph_MapColoring::OriginalTypeRef orig)
Generates the Epetra_MapColoring object from an input Epetra_CrsGraph.
int NumMyElements() const
int SetSeed(unsigned int Seed_in)
const Epetra_BlockMap & RowMap() const
const Epetra_Comm & Comm() const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int NumColors() const
int NumMyIndices(int Row) const
int MaxNumIndices() const
Map Coloring of independent columns in a Graph.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)