Ifpack_METISReordering.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #ifdef HAVE_IFPACK_TEUCHOS
00003 #include "Ifpack_Reordering.h"
00004 #include "Ifpack_METISReordering.h"
00005 #include "Ifpack_Graph.h"
00006 #include "Ifpack_Graph_Epetra_CrsGraph.h"
00007 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00008 #include "Epetra_Comm.h"
00009 #include "Epetra_MultiVector.h"
00010 #include "Epetra_CrsGraph.h"
00011 #include "Epetra_Map.h"
00012 #include "Teuchos_ParameterList.hpp"
00013 
00014 typedef int idxtype;
00015 #ifdef HAVE_IFPACK_METIS
00016 extern "C" {
00017   void METIS_NodeND(int *n, idxtype *xadj, idxtype *adjncy, 
00018             int *numflag, int *options, int *perm, int *iperm);
00019 }
00020 #endif
00021 
00022 //==============================================================================
00023 Ifpack_METISReordering::Ifpack_METISReordering() :
00024   UseSymmetricGraph_(false),
00025   NumMyRows_(0),
00026   IsComputed_(false)
00027 { }
00028 
00029 //==============================================================================
00030 // Mainly copied from Ifpack_METISPartitioner.cpp
00031 //
00032 // NOTE:
00033 // - matrix is supposed to be localized, and passes through the
00034 // singleton filter. This means that I do not have to look
00035 // for Dirichlet nodes (singletons). Also, all rows and columns are 
00036 // local.
00037 int Ifpack_METISReordering::Compute(const Ifpack_Graph& Graph)
00038 {
00039 
00040   NumMyRows_ = Graph.NumMyRows();
00041   Reorder_.resize(NumMyRows_);
00042   InvReorder_.resize(NumMyRows_);
00043 
00044   int ierr;
00045 
00046   Epetra_CrsGraph* SymGraph = 0;
00047   Epetra_Map* SymMap = 0;
00048   Ifpack_Graph_Epetra_CrsGraph* SymIFPACKGraph = 0;
00049   Ifpack_Graph* IFPACKGraph = (Ifpack_Graph*)&Graph;
00050 
00051   int Length = 2 * Graph.MaxMyNumEntries();
00052   int NumIndices;
00053   vector<int> Indices;
00054   Indices.resize(Length);
00055 
00056   vector<int> options;
00057   options.resize(8);
00058   options[0] = 0; // default values
00059 
00060 #ifdef HAVE_IFPACK_METIS  
00061   int numflag = 0; // C style
00062 #endif
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           // insert A(i,j), then A(j,i)
00083       SymGraph->InsertGlobalIndices(i,1,&jj);
00084       SymGraph->InsertGlobalIndices(jj,1,&i);
00085     }
00086       }      
00087     }
00088     IFPACK_CHK_ERR(SymGraph->OptimizeStorage());
00089     IFPACK_CHK_ERR(SymGraph->FillComplete());
00090     SymIFPACKGraph = new Ifpack_Graph_Epetra_CrsGraph(SymGraph);
00091     IFPACKGraph = SymIFPACKGraph;
00092   }
00093 
00094   // convert to METIS format
00095   vector<idxtype> xadj;
00096   xadj.resize(NumMyRows_ + 1);
00097 
00098   vector<idxtype> adjncy;
00099   adjncy.resize(Graph.NumMyNonzeros());
00100    
00101   int count = 0; 
00102   int count2 = 0; 
00103   xadj[0] = 0;
00104   
00105   for (int i = 0; i < NumMyRows_ ; ++i) {
00106 
00107     xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
00108 
00109     ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
00110     IFPACK_CHK_ERR(ierr);
00111 
00112     for (int j = 0 ; j < NumIndices ; ++j) {
00113       int jj = Indices[j];
00114       if (jj != i) {
00115     adjncy[count++] = jj;
00116     xadj[count2+1]++;
00117       }
00118     }
00119     count2++;
00120   }
00121 
00122 #ifdef HAVE_IFPACK_METIS
00123   // vectors from METIS. The second last is `perm', the last is `iperm'.
00124   // They store the fill-reducing permutation and inverse-permutation.
00125   // Let A be the original matrix and A' the permuted matrix. The
00126   // arrays perm and iperm are defined as follows. Row (column) i of A'
00127   // if the perm[i] row (col) of A, and row (column) i of A is the
00128   // iperm[i] row (column) of A'. The numbering starts from 0 in our case.
00129   METIS_NodeND(&NumMyRows_, &xadj[0], &adjncy[0],
00130            &numflag, &options[0],
00131            &InvReorder_[0], &Reorder_[0]);
00132 #else
00133   cerr << "Please configure with --enable-ifpack-metis" << endl;
00134   cerr << "to use METIS Reordering." << endl;
00135   exit(EXIT_FAILURE);
00136 #endif
00137       
00138   if (SymGraph)
00139     delete SymGraph;
00140   if (SymMap)
00141     delete SymMap;
00142   if (SymIFPACKGraph)
00143     delete SymIFPACKGraph;
00144 
00145   return(0);
00146 } 
00147 
00148 //==============================================================================
00149 int Ifpack_METISReordering::Compute(const Epetra_RowMatrix& Matrix)
00150 {
00151   Ifpack_Graph_Epetra_RowMatrix Graph(&Matrix);
00152 
00153   IFPACK_CHK_ERR(Compute(Graph));
00154 
00155   return(0);
00156 }
00157 
00158 //==============================================================================
00159 int Ifpack_METISReordering::Reorder(const int i) const
00160 {
00161 #ifdef IFPACK_ABC
00162   if (!IsComputed())
00163     IFPACK_CHK_ERR(-1);
00164   if ((i < 0) || (i >= NumMyRows_))
00165     IFPACK_CHK_ERR(-1);
00166 #endif
00167 
00168   return(Reorder_[i]);
00169 }
00170 
00171 //==============================================================================
00172 int Ifpack_METISReordering::InvReorder(const int i) const
00173 {
00174 #ifdef IFPACK_ABC
00175   if (!IsComputed())
00176     IFPACK_CHK_ERR(-1);
00177   if ((i < 0) || (i >= NumMyRows_))
00178     IFPACK_CHK_ERR(-1);
00179 #endif
00180 
00181   return(InvReorder_[i]);
00182 }
00183 //==============================================================================
00184 int Ifpack_METISReordering::P(const Epetra_MultiVector& Xorig,
00185                 Epetra_MultiVector& X) const
00186 {  
00187   int NumVectors = X.NumVectors();
00188 
00189   for (int j = 0 ; j < NumVectors ; ++j) {
00190     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00191       int np = Reorder_[i];
00192       X[j][np] = Xorig[j][i];
00193     }
00194   }
00195 
00196   return(0);
00197 }
00198 
00199 //==============================================================================
00200 int Ifpack_METISReordering::Pinv(const Epetra_MultiVector& Xorig,
00201                  Epetra_MultiVector& X) const
00202 {
00203   int NumVectors = X.NumVectors();
00204 
00205   for (int j = 0 ; j < NumVectors ; ++j) {
00206     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00207       int np = Reorder_[i];
00208       X[j][i] = Xorig[j][np];
00209     }
00210   }
00211 
00212   return(0);
00213 }
00214 
00215 //==============================================================================
00216 ostream& Ifpack_METISReordering::Print(std::ostream& os) const
00217 {
00218   os << "*** Ifpack_METISReordering" << endl << endl;
00219   if (!IsComputed())
00220     os << "*** Reordering not yet computed." << endl;
00221   
00222   os << "*** Number of local rows = " << NumMyRows_ << endl;
00223   os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
00224   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00225     os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
00226   }
00227    
00228   return(os);
00229 }
00230 
00231 #endif

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