Ifpack_METISReordering.cpp

Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Ifpack_ConfigDefs.h"
00031 #include "Ifpack_Reordering.h"
00032 #include "Ifpack_METISReordering.h"
00033 #include "Ifpack_Graph.h"
00034 #include "Ifpack_Graph_Epetra_CrsGraph.h"
00035 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00036 #include "Epetra_Comm.h"
00037 #include "Epetra_MultiVector.h"
00038 #include "Epetra_CrsGraph.h"
00039 #include "Epetra_Map.h"
00040 #include "Teuchos_ParameterList.hpp"
00041 
00042 typedef int idxtype;
00043 #ifdef HAVE_IFPACK_METIS
00044 extern "C" {
00045   void METIS_NodeND(int *n, idxtype *xadj, idxtype *adjncy, 
00046         int *numflag, int *options, int *perm, int *iperm);
00047 }
00048 #endif
00049 
00050 //==============================================================================
00051 Ifpack_METISReordering::Ifpack_METISReordering() :
00052   UseSymmetricGraph_(false),
00053   NumMyRows_(0),
00054   IsComputed_(false)
00055 {}
00056 
00057 //==============================================================================
00058 // Mainly copied from Ifpack_METISPartitioner.cpp
00059 //
00060 // NOTE:
00061 // - matrix is supposed to be localized, and passes through the
00062 // singleton filter. This means that I do not have to look
00063 // for Dirichlet nodes (singletons). Also, all rows and columns are 
00064 // local.
00065 int Ifpack_METISReordering::Compute(const Ifpack_Graph& Graph)
00066 {
00067 
00068   NumMyRows_ = Graph.NumMyRows();
00069   Reorder_.resize(NumMyRows_);
00070   InvReorder_.resize(NumMyRows_);
00071 
00072   int ierr;
00073 
00074   Teuchos::RefCountPtr<Epetra_CrsGraph> SymGraph;
00075   Teuchos::RefCountPtr<Epetra_Map> SymMap;
00076   Teuchos::RefCountPtr<Ifpack_Graph_Epetra_CrsGraph> SymIFPACKGraph;
00077   Teuchos::RefCountPtr<Ifpack_Graph> IFPACKGraph = Teuchos::rcp( (Ifpack_Graph*)&Graph, false );
00078 
00079   int Length = 2 * Graph.MaxMyNumEntries();
00080   int NumIndices;
00081   vector<int> Indices;
00082   Indices.resize(Length);
00083 
00084   vector<int> options;
00085   options.resize(8);
00086   options[0] = 0; // default values
00087 
00088 #ifdef HAVE_IFPACK_METIS  
00089   int numflag = 0; // C style
00090 #endif
00091 
00092   if (UseSymmetricGraph_) {
00093 
00094     // need to build a symmetric graph. 
00095     // I do this in two stages:
00096     // 1.- construct an Epetra_CrsMatrix, symmetric
00097     // 2.- convert the Epetra_CrsMatrix into METIS format
00098     SymMap = Teuchos::rcp( new Epetra_Map(NumMyRows_,0,Graph.Comm()) );
00099     SymGraph = Teuchos::rcp( new Epetra_CrsGraph(Copy,*SymMap,0) );
00100 
00101     for (int i = 0; i < NumMyRows_ ; ++i) {
00102 
00103       ierr = Graph.ExtractMyRowCopy(i, Length, NumIndices, 
00104               &Indices[0]);
00105       IFPACK_CHK_ERR(ierr);
00106 
00107       for (int j = 0 ; j < NumIndices ; ++j) {
00108   int jj = Indices[j];
00109   if (jj != i) {
00110           // insert A(i,j), then A(j,i)
00111     SymGraph->InsertGlobalIndices(i,1,&jj);
00112     SymGraph->InsertGlobalIndices(jj,1,&i);
00113   }
00114       }      
00115     }
00116     IFPACK_CHK_ERR(SymGraph->OptimizeStorage());
00117     IFPACK_CHK_ERR(SymGraph->FillComplete());
00118     SymIFPACKGraph = Teuchos::rcp( new Ifpack_Graph_Epetra_CrsGraph(SymGraph) );
00119     IFPACKGraph = SymIFPACKGraph;
00120   }
00121 
00122   // convert to METIS format
00123   vector<idxtype> xadj;
00124   xadj.resize(NumMyRows_ + 1);
00125 
00126   vector<idxtype> adjncy;
00127   adjncy.resize(Graph.NumMyNonzeros());
00128    
00129   int count = 0; 
00130   int count2 = 0; 
00131   xadj[0] = 0;
00132   
00133   for (int i = 0; i < NumMyRows_ ; ++i) {
00134 
00135     xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
00136 
00137     ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
00138     IFPACK_CHK_ERR(ierr);
00139 
00140     for (int j = 0 ; j < NumIndices ; ++j) {
00141       int jj = Indices[j];
00142       if (jj != i) {
00143   adjncy[count++] = jj;
00144   xadj[count2+1]++;
00145       }
00146     }
00147     count2++;
00148   }
00149 
00150 #ifdef HAVE_IFPACK_METIS
00151   // vectors from METIS. The second last is `perm', the last is `iperm'.
00152   // They store the fill-reducing permutation and inverse-permutation.
00153   // Let A be the original matrix and A' the permuted matrix. The
00154   // arrays perm and iperm are defined as follows. Row (column) i of A'
00155   // if the perm[i] row (col) of A, and row (column) i of A is the
00156   // iperm[i] row (column) of A'. The numbering starts from 0 in our case.
00157   METIS_NodeND(&NumMyRows_, &xadj[0], &adjncy[0],
00158          &numflag, &options[0],
00159          &InvReorder_[0], &Reorder_[0]);
00160 #else
00161   cerr << "Please configure with --enable-ifpack-metis" << endl;
00162   cerr << "to use METIS Reordering." << endl;
00163   exit(EXIT_FAILURE);
00164 #endif
00165       
00166   return(0);
00167 } 
00168 
00169 //==============================================================================
00170 int Ifpack_METISReordering::Compute(const Epetra_RowMatrix& Matrix)
00171 {
00172   Ifpack_Graph_Epetra_RowMatrix Graph(Teuchos::rcp(&Matrix, false));
00173 
00174   IFPACK_CHK_ERR(Compute(Graph));
00175 
00176   return(0);
00177 }
00178 
00179 //==============================================================================
00180 int Ifpack_METISReordering::Reorder(const int i) const
00181 {
00182 #ifdef IFPACK_ABC
00183   if (!IsComputed())
00184     IFPACK_CHK_ERR(-1);
00185   if ((i < 0) || (i >= NumMyRows_))
00186     IFPACK_CHK_ERR(-1);
00187 #endif
00188 
00189   return(Reorder_[i]);
00190 }
00191 
00192 //==============================================================================
00193 int Ifpack_METISReordering::InvReorder(const int i) const
00194 {
00195 #ifdef IFPACK_ABC
00196   if (!IsComputed())
00197     IFPACK_CHK_ERR(-1);
00198   if ((i < 0) || (i >= NumMyRows_))
00199     IFPACK_CHK_ERR(-1);
00200 #endif
00201 
00202   return(InvReorder_[i]);
00203 }
00204 //==============================================================================
00205 int Ifpack_METISReordering::P(const Epetra_MultiVector& Xorig,
00206           Epetra_MultiVector& X) const
00207 {  
00208   int NumVectors = X.NumVectors();
00209 
00210   for (int j = 0 ; j < NumVectors ; ++j) {
00211     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00212       int np = Reorder_[i];
00213       X[j][np] = Xorig[j][i];
00214     }
00215   }
00216 
00217   return(0);
00218 }
00219 
00220 //==============================================================================
00221 int Ifpack_METISReordering::Pinv(const Epetra_MultiVector& Xorig,
00222          Epetra_MultiVector& X) const
00223 {
00224   int NumVectors = X.NumVectors();
00225 
00226   for (int j = 0 ; j < NumVectors ; ++j) {
00227     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00228       int np = Reorder_[i];
00229       X[j][i] = Xorig[j][np];
00230     }
00231   }
00232 
00233   return(0);
00234 }
00235 
00236 //==============================================================================
00237 ostream& Ifpack_METISReordering::Print(std::ostream& os) const
00238 {
00239   os << "*** Ifpack_METISReordering" << endl << endl;
00240   if (!IsComputed())
00241     os << "*** Reordering not yet computed." << endl;
00242   
00243   os << "*** Number of local rows = " << NumMyRows_ << endl;
00244   os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
00245   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00246     os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
00247   }
00248    
00249   return(os);
00250 }
00251 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:35 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3