IFPACK Development
Ifpack_METISReordering.cpp
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #include "Ifpack_ConfigDefs.h"
00044 #include "Ifpack_Reordering.h"
00045 #include "Ifpack_METISReordering.h"
00046 #include "Ifpack_Graph.h"
00047 #include "Ifpack_Graph_Epetra_CrsGraph.h"
00048 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00049 #include "Epetra_Comm.h"
00050 #include "Epetra_MultiVector.h"
00051 #include "Epetra_CrsGraph.h"
00052 #include "Epetra_Map.h"
00053 #include "Teuchos_ParameterList.hpp"
00054 
00055 typedef int idxtype;
00056 #ifdef HAVE_IFPACK_METIS
00057 extern "C" {
00058   void METIS_NodeND(int *n, idxtype *xadj, idxtype *adjncy, 
00059             int *numflag, int *options, int *perm, int *iperm);
00060 }
00061 #endif
00062 
00063 //==============================================================================
00064 Ifpack_METISReordering::Ifpack_METISReordering() :
00065   UseSymmetricGraph_(false),
00066   NumMyRows_(0),
00067   IsComputed_(false)
00068 {}
00069 
00070 //==============================================================================
00071 // Mainly copied from Ifpack_METISPartitioner.cpp
00072 //
00073 // NOTE:
00074 // - matrix is supposed to be localized, and passes through the
00075 // singleton filter. This means that I do not have to look
00076 // for Dirichlet nodes (singletons). Also, all rows and columns are 
00077 // local.
00078 int Ifpack_METISReordering::Compute(const Ifpack_Graph& Graph)
00079 {
00080 
00081   NumMyRows_ = Graph.NumMyRows();
00082   Reorder_.resize(NumMyRows_);
00083   InvReorder_.resize(NumMyRows_);
00084 
00085   int ierr;
00086 
00087   Teuchos::RefCountPtr<Epetra_CrsGraph> SymGraph;
00088   Teuchos::RefCountPtr<Epetra_Map> SymMap;
00089   Teuchos::RefCountPtr<Ifpack_Graph_Epetra_CrsGraph> SymIFPACKGraph;
00090   Teuchos::RefCountPtr<Ifpack_Graph> IFPACKGraph = Teuchos::rcp( (Ifpack_Graph*)&Graph, false );
00091 
00092   int Length = 2 * Graph.MaxMyNumEntries();
00093   int NumIndices;
00094   std::vector<int> Indices;
00095   Indices.resize(Length);
00096 
00097   std::vector<int> options;
00098   options.resize(8);
00099   options[0] = 0; // default values
00100 
00101 #ifdef HAVE_IFPACK_METIS  
00102   int numflag = 0; // C style
00103 #endif
00104 
00105   if (UseSymmetricGraph_) {
00106 
00107 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00108     // need to build a symmetric graph. 
00109     // I do this in two stages:
00110     // 1.- construct an Epetra_CrsMatrix, symmetric
00111     // 2.- convert the Epetra_CrsMatrix into METIS format
00112     SymMap = Teuchos::rcp( new Epetra_Map(NumMyRows_,0,Graph.Comm()) );
00113     SymGraph = Teuchos::rcp( new Epetra_CrsGraph(Copy,*SymMap,0) );
00114 #endif
00115 
00116 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00117     if(SymGraph->RowMap().GlobalIndicesInt()) {
00118     for (int i = 0; i < NumMyRows_ ; ++i) {
00119 
00120       ierr = Graph.ExtractMyRowCopy(i, Length, NumIndices, 
00121                       &Indices[0]);
00122       IFPACK_CHK_ERR(ierr);
00123 
00124       for (int j = 0 ; j < NumIndices ; ++j) {
00125     int jj = Indices[j];
00126     if (jj != i) {
00127           // insert A(i,j), then A(j,i)
00128       SymGraph->InsertGlobalIndices(i,1,&jj);
00129       SymGraph->InsertGlobalIndices(jj,1,&i);
00130     }
00131       }      
00132     }
00133     }
00134     else
00135 #endif
00136 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00137     if(SymGraph->RowMap().GlobalIndicesLongLong()) {
00138     for (int i = 0; i < NumMyRows_ ; ++i) {
00139       long long i_LL = i;
00140 
00141       ierr = Graph.ExtractMyRowCopy(i, Length, NumIndices, 
00142                       &Indices[0]);
00143       IFPACK_CHK_ERR(ierr);
00144 
00145       for (int j = 0 ; j < NumIndices ; ++j) {
00146     long long jj = Indices[j];
00147     if (jj != i) {
00148           // insert A(i,j), then A(j,i)
00149       SymGraph->InsertGlobalIndices(i_LL,1,&jj);
00150       SymGraph->InsertGlobalIndices(jj,1,&i_LL);
00151     }
00152       }      
00153     }
00154     }
00155     else
00156 #endif
00157     throw "Ifpack_METISReordering::Compute: GlobalIndices type unknown";
00158 
00159     IFPACK_CHK_ERR(SymGraph->OptimizeStorage());
00160     IFPACK_CHK_ERR(SymGraph->FillComplete());
00161     SymIFPACKGraph = Teuchos::rcp( new Ifpack_Graph_Epetra_CrsGraph(SymGraph) );
00162     IFPACKGraph = SymIFPACKGraph;
00163   }
00164 
00165   // convert to METIS format
00166   std::vector<idxtype> xadj;
00167   xadj.resize(NumMyRows_ + 1);
00168 
00169   std::vector<idxtype> adjncy;
00170   adjncy.resize(Graph.NumMyNonzeros());
00171    
00172   int count = 0; 
00173   int count2 = 0; 
00174   xadj[0] = 0;
00175   
00176   for (int i = 0; i < NumMyRows_ ; ++i) {
00177 
00178     xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
00179 
00180     ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
00181     IFPACK_CHK_ERR(ierr);
00182 
00183     for (int j = 0 ; j < NumIndices ; ++j) {
00184       int jj = Indices[j];
00185       if (jj != i) {
00186     adjncy[count++] = jj;
00187     xadj[count2+1]++;
00188       }
00189     }
00190     count2++;
00191   }
00192 
00193 #ifdef HAVE_IFPACK_METIS
00194   // vectors from METIS. The second last is `perm', the last is `iperm'.
00195   // They store the fill-reducing permutation and inverse-permutation.
00196   // Let A be the original matrix and A' the permuted matrix. The
00197   // arrays perm and iperm are defined as follows. Row (column) i of A'
00198   // if the perm[i] row (col) of A, and row (column) i of A is the
00199   // iperm[i] row (column) of A'. The numbering starts from 0 in our case.
00200   METIS_NodeND(&NumMyRows_, &xadj[0], &adjncy[0],
00201            &numflag, &options[0],
00202            &InvReorder_[0], &Reorder_[0]);
00203 #else
00204   cerr << "Please configure with --enable-ifpack-metis" << endl;
00205   cerr << "to use METIS Reordering." << endl;
00206   exit(EXIT_FAILURE);
00207 #endif
00208       
00209   return(0);
00210 } 
00211 
00212 //==============================================================================
00213 int Ifpack_METISReordering::Compute(const Epetra_RowMatrix& Matrix)
00214 {
00215   Ifpack_Graph_Epetra_RowMatrix Graph(Teuchos::rcp(&Matrix, false));
00216 
00217   IFPACK_CHK_ERR(Compute(Graph));
00218 
00219   return(0);
00220 }
00221 
00222 //==============================================================================
00223 int Ifpack_METISReordering::Reorder(const int i) const
00224 {
00225 #ifdef IFPACK_ABC
00226   if (!IsComputed())
00227     IFPACK_CHK_ERR(-1);
00228   if ((i < 0) || (i >= NumMyRows_))
00229     IFPACK_CHK_ERR(-1);
00230 #endif
00231 
00232   return(Reorder_[i]);
00233 }
00234 
00235 //==============================================================================
00236 int Ifpack_METISReordering::InvReorder(const int i) const
00237 {
00238 #ifdef IFPACK_ABC
00239   if (!IsComputed())
00240     IFPACK_CHK_ERR(-1);
00241   if ((i < 0) || (i >= NumMyRows_))
00242     IFPACK_CHK_ERR(-1);
00243 #endif
00244 
00245   return(InvReorder_[i]);
00246 }
00247 //==============================================================================
00248 int Ifpack_METISReordering::P(const Epetra_MultiVector& Xorig,
00249                 Epetra_MultiVector& X) const
00250 {  
00251   int NumVectors = X.NumVectors();
00252 
00253   for (int j = 0 ; j < NumVectors ; ++j) {
00254     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00255       int np = Reorder_[i];
00256       X[j][np] = Xorig[j][i];
00257     }
00258   }
00259 
00260   return(0);
00261 }
00262 
00263 //==============================================================================
00264 int Ifpack_METISReordering::Pinv(const Epetra_MultiVector& Xorig,
00265                  Epetra_MultiVector& X) const
00266 {
00267   int NumVectors = X.NumVectors();
00268 
00269   for (int j = 0 ; j < NumVectors ; ++j) {
00270     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00271       int np = Reorder_[i];
00272       X[j][i] = Xorig[j][np];
00273     }
00274   }
00275 
00276   return(0);
00277 }
00278 
00279 //==============================================================================
00280 ostream& Ifpack_METISReordering::Print(std::ostream& os) const
00281 {
00282   os << "*** Ifpack_METISReordering" << endl << endl;
00283   if (!IsComputed())
00284     os << "*** Reordering not yet computed." << endl;
00285   
00286   os << "*** Number of local rows = " << NumMyRows_ << endl;
00287   os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
00288   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00289     os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
00290   }
00291    
00292   return(os);
00293 }
00294 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends