Ifpack_METISPartitioner.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #ifdef HAVE_IFPACK_TEUCHOS
00003 #include "Ifpack_Partitioner.h"
00004 #include "Ifpack_OverlappingPartitioner.h"
00005 #include "Ifpack_METISPartitioner.h"
00006 #include "Ifpack_Graph.h"
00007 #include "Ifpack_Graph_Epetra_CrsGraph.h"
00008 #include "Epetra_Comm.h"
00009 #include "Epetra_CrsGraph.h"
00010 #include "Epetra_Map.h"
00011 #include "Teuchos_ParameterList.hpp"
00012 
00013 // may need to change this for wierd installations
00014 typedef int idxtype;
00015 #ifdef HAVE_IFPACK_METIS
00016 extern "C" {
00017   void METIS_EstimateMemory(int *, idxtype *, idxtype *, int *, int *, int *);
00018   void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, 
00019                idxtype *, int *, int *, int *, int *, int *,
00020                idxtype *);
00021   void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, 
00022                 idxtype *, idxtype *, int *, int *, int *, 
00023                 int *, int *, idxtype *);
00024 
00025 }
00026 #endif
00027 
00028 //==============================================================================
00029 // NOTE:
00030 // - matrix is supposed to be localized, and passes through the
00031 // singleton filter. This means that I do not have to look
00032 // for Dirichlet nodes (singletons). Also, all rows and columns are 
00033 // local.
00034 int Ifpack_METISPartitioner::ComputePartitions()
00035 {
00036 
00037   int ierr;
00038 #ifdef HAVE_IFPACK_METIS
00039   int nbytes = 0;
00040   int edgecut;
00041 #endif
00042 
00043   Epetra_CrsGraph* SymGraph = 0;
00044   Epetra_Map* SymMap = 0;
00045   Ifpack_Graph_Epetra_CrsGraph* SymIFPACKGraph = 0;
00046   Ifpack_Graph* IFPACKGraph = (Ifpack_Graph*)Graph_;
00047 
00048   int Length = 2 * MaxNumEntries();
00049   int NumIndices;
00050   vector<int> Indices;
00051   Indices.resize(Length);
00052 
00053   /* construct the CSR graph information of the LOCAL matrix
00054      using the get_row function */
00055 
00056   vector<idxtype> wgtflag;
00057   wgtflag.resize(4);
00058 
00059   vector<int> options;
00060   options.resize(4);
00061   
00062   int numflag;
00063 
00064   if (UseSymmetricGraph_) {
00065 
00066     // need to build a symmetric graph. 
00067     // I do this in two stages:
00068     // 1.- construct an Epetra_CrsMatrix, symmetric
00069     // 2.- convert the Epetra_CrsMatrix into METIS format
00070     SymMap = new Epetra_Map(NumMyRows(),0,Graph_->Comm());
00071     SymGraph = new Epetra_CrsGraph(Copy,*SymMap,0);
00072 
00073     for (int i = 0; i < NumMyRows() ; ++i) {
00074 
00075       ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, 
00076                       &Indices[0]);
00077       IFPACK_CHK_ERR(ierr);
00078 
00079       for (int j = 0 ; j < NumIndices ; ++j) {
00080     int jj = Indices[j];
00081     if (jj != i) {
00082       SymGraph->InsertGlobalIndices(i,1,&jj);
00083       SymGraph->InsertGlobalIndices(jj,1,&i);
00084     }
00085       }      
00086     }
00087     IFPACK_CHK_ERR(SymGraph->FillComplete());
00088     SymIFPACKGraph = new Ifpack_Graph_Epetra_CrsGraph(SymGraph);
00089     IFPACKGraph = SymIFPACKGraph;
00090   }
00091 
00092   // now work on IFPACKGraph, that can be the symmetric or
00093   // the non-symmetric one
00094 
00095   /* set parameters */
00096    
00097   wgtflag[0] = 0;    /* no weights */
00098   numflag    = 0;    /* C style */
00099   options[0] = 0;    /* default options */
00100    
00101   vector<idxtype> xadj;
00102   xadj.resize(NumMyRows() + 1);
00103 
00104   vector<idxtype> adjncy;
00105   adjncy.resize(NumMyNonzeros());
00106    
00107   int count = 0; 
00108   int count2 = 0; 
00109   xadj[0] = 0;
00110   
00111   for (int i = 0; i < NumMyRows() ; ++i) {
00112 
00113     xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
00114 
00115     ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
00116     IFPACK_CHK_ERR(ierr);
00117 
00118     for (int j = 0 ; j < NumIndices ; ++j) {
00119       int jj = Indices[j];
00120       if (jj != i) {
00121     adjncy[count++] = jj;
00122     xadj[count2+1]++;
00123       }
00124     }
00125     count2++;
00126   }
00127 
00128   vector<idxtype> NodesInSubgraph;
00129   NodesInSubgraph.resize(NumLocalParts_);
00130 
00131   // some cases can be handled separately
00132   
00133   int ok;
00134 
00135   if (NumLocalParts() == 1) {
00136 
00137     for (int i = 0 ; i < NumMyRows() ; ++i) 
00138       Partition_[i] = 0;
00139     
00140   } else if (NumLocalParts() == NumMyRows()) {
00141 
00142     for (int i = 0 ; i < NumMyRows() ; ++i) 
00143       Partition_[i] = i;
00144   
00145   } else {
00146 
00147     ok = 0;
00148 
00149     // sometimes METIS creates less partitions than specified.
00150     // ok will check this problem, and recall metis, asking
00151     // for NumLocalParts_/2 partitions
00152     while (ok == 0) {
00153       
00154       for (int i = 0 ; i < NumMyRows() ; ++i) 
00155     Partition_[i] = -1;
00156     
00157 #ifdef HAVE_IFPACK_METIS
00158       int j = NumMyRows();
00159       if (NumLocalParts_ < 8) {
00160 
00161     int i = 1; /* optype in the METIS manual */
00162     numflag = 0;
00163     METIS_EstimateMemory(&j, &xadj[0], &adjncy[0], 
00164                  &numflag, &i, &nbytes );
00165     
00166     METIS_PartGraphRecursive(&j, &xadj[0], &adjncy[0],
00167                  NULL, NULL,
00168                  &wgtflag[0], &numflag, &NumLocalParts_, 
00169                  &options[0], &edgecut, &Partition_[0]);
00170       } else {
00171 
00172     numflag = 0;
00173     
00174     METIS_PartGraphKway (&j, &xadj[0], &adjncy[0], 
00175                  NULL, 
00176                  NULL, &wgtflag[0], &numflag, 
00177                  &NumLocalParts_, &options[0],
00178                  &edgecut, &Partition_[0]);
00179       }
00180 #else
00181       numflag = numflag * 2; // avoid warning for unused variable
00182       if (Graph_->Comm().MyPID() == 0) {
00183     cerr << "METIS was not linked; now I put all" << endl;
00184     cerr << "the local nodes in the same partition." << endl;
00185       }
00186       for (int i = 0 ; i < NumMyRows() ; ++i) 
00187     Partition_[i] = 0;
00188       NumLocalParts_ = 1;
00189 #endif
00190       
00191       ok = 1;
00192       
00193       for (int i = 0 ; i < NumLocalParts() ; ++i) 
00194     NodesInSubgraph[i] = 0;
00195 
00196       for (int i = 0 ; i < NumMyRows() ; ++i) {
00197     int j = Partition_[i];
00198     if ((j < 0) || (j>= NumLocalParts())) {
00199       ok = 0;
00200       break;
00201     } 
00202     else NodesInSubgraph[j]++;
00203       }
00204       
00205       for (int i = 0 ; i < NumLocalParts() ; ++i) {
00206     if( NodesInSubgraph[i] == 0 ) {
00207       ok = 0;
00208       break;
00209     }
00210       }
00211       
00212       if (ok == 0) {
00213     cerr << "Specified number of subgraphs ("
00214          << NumLocalParts_ << ") generates empty subgraphs." << endl;
00215     cerr << "Now I recall METIS with NumLocalParts_ = "
00216          << NumLocalParts_ / 2 << "..." << endl;
00217     NumLocalParts_ = NumLocalParts_/2;
00218       }
00219       
00220       if (NumLocalParts() == 0) {
00221     IFPACK_CHK_ERR(-10); // something went wrong
00222       }
00223       
00224       if (NumLocalParts() == 1) {
00225     for (int i = 0 ; i < NumMyRows() ; ++i) 
00226       Partition_[i] = 0;
00227     ok = 1;
00228       }
00229       
00230     } /* while( ok == 0 ) */
00231   
00232   } /* if( NumLocalParts_ == 1 ) */
00233 
00234   if (SymGraph)
00235     delete SymGraph;
00236   if (SymMap)
00237     delete SymMap;
00238   if (SymIFPACKGraph)
00239     delete SymIFPACKGraph;
00240 
00241   return(0);
00242 } 
00243 #endif

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1