IFPACK Development
Ifpack_METISPartitioner.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_Partitioner.h"
00045 #include "Ifpack_OverlappingPartitioner.h"
00046 #include "Ifpack_METISPartitioner.h"
00047 #include "Ifpack_Graph.h"
00048 #include "Ifpack_Graph_Epetra_CrsGraph.h"
00049 #include "Epetra_Comm.h"
00050 #include "Epetra_CrsGraph.h"
00051 #include "Epetra_Map.h"
00052 #include "Teuchos_ParameterList.hpp"
00053 #include "Teuchos_RefCountPtr.hpp"
00054 
00055 // may need to change this for wierd installations
00056 typedef int idxtype;
00057 #ifdef HAVE_IFPACK_METIS
00058 extern "C" {
00059   void METIS_EstimateMemory(int *, idxtype *, idxtype *, int *, int *, int *);
00060   void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, 
00061                idxtype *, int *, int *, int *, int *, int *,
00062                idxtype *);
00063   void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, 
00064                 idxtype *, idxtype *, int *, int *, int *, 
00065                 int *, int *, idxtype *);
00066 
00067 }
00068 #endif
00069 
00070 //==============================================================================
00071 // NOTE:
00072 // - matrix is supposed to be localized, and passes through the
00073 // singleton filter. This means that I do not have to look
00074 // for Dirichlet nodes (singletons). Also, all rows and columns are 
00075 // local.
00076 int Ifpack_METISPartitioner::ComputePartitions()
00077 {
00078 
00079   int ierr;
00080 #ifdef HAVE_IFPACK_METIS
00081   int nbytes = 0;
00082   int edgecut;
00083 #endif
00084 
00085   Teuchos::RefCountPtr<Epetra_CrsGraph> SymGraph ;
00086   Teuchos::RefCountPtr<Epetra_Map> SymMap;
00087   Teuchos::RefCountPtr<Ifpack_Graph_Epetra_CrsGraph> SymIFPACKGraph;
00088   Teuchos::RefCountPtr<Ifpack_Graph> IFPACKGraph = Teuchos::rcp( (Ifpack_Graph*)Graph_, false );
00089 
00090   int Length = 2 * MaxNumEntries();
00091   int NumIndices;
00092   std::vector<int> Indices;
00093   Indices.resize(Length);
00094 
00095   /* construct the CSR graph information of the LOCAL matrix
00096      using the get_row function */
00097 
00098   std::vector<idxtype> wgtflag;
00099   wgtflag.resize(4);
00100 
00101   std::vector<int> options;
00102   options.resize(4);
00103   
00104   int numflag;
00105 
00106   if (UseSymmetricGraph_) {
00107 
00108 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00109     // need to build a symmetric graph. 
00110     // I do this in two stages:
00111     // 1.- construct an Epetra_CrsMatrix, symmetric
00112     // 2.- convert the Epetra_CrsMatrix into METIS format
00113     SymMap = Teuchos::rcp( new Epetra_Map(NumMyRows(),0,Graph_->Comm()) );
00114     SymGraph = Teuchos::rcp( new Epetra_CrsGraph(Copy,*SymMap,0) );
00115 #endif
00116 
00117 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00118       if(SymGraph->RowMap().GlobalIndicesInt()) {
00119         for (int i = 0; i < NumMyRows() ; ++i) {
00120 
00121           ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, &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               SymGraph->InsertGlobalIndices(i,1,&jj);
00128               SymGraph->InsertGlobalIndices(jj,1,&i);
00129             }
00130           }
00131         }
00132       }
00133       else
00134 #endif
00135 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00136       if(SymGraph->RowMap().GlobalIndicesLongLong()) {
00137         for (int i = 0; i < NumMyRows() ; ++i) {
00138           long long i_LL = i;
00139 
00140           ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
00141           IFPACK_CHK_ERR(ierr);
00142 
00143           for (int j = 0 ; j < NumIndices ; ++j) {
00144             long long jj = Indices[j];
00145             if (jj != i_LL) {
00146               SymGraph->InsertGlobalIndices(i_LL,1,&jj);
00147               SymGraph->InsertGlobalIndices(jj,1,&i_LL);
00148             }
00149           }
00150         }
00151       }
00152       else
00153 #endif
00154         throw "Ifpack_METISPartitioner::ComputePartitions: GlobalIndices type unknown";
00155 
00156     IFPACK_CHK_ERR(SymGraph->FillComplete());
00157     SymIFPACKGraph = Teuchos::rcp( new Ifpack_Graph_Epetra_CrsGraph(SymGraph) );
00158     IFPACKGraph = SymIFPACKGraph;
00159   }
00160 
00161   // now work on IFPACKGraph, that can be the symmetric or
00162   // the non-symmetric one
00163 
00164   /* set parameters */
00165    
00166   wgtflag[0] = 0;    /* no weights */
00167   numflag    = 0;    /* C style */
00168   options[0] = 0;    /* default options */
00169    
00170   std::vector<idxtype> xadj;
00171   xadj.resize(NumMyRows() + 1);
00172 
00173   std::vector<idxtype> adjncy;
00174   adjncy.resize(NumMyNonzeros());
00175    
00176   int count = 0; 
00177   int count2 = 0; 
00178   xadj[0] = 0;
00179   
00180   for (int i = 0; i < NumMyRows() ; ++i) {
00181 
00182     xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
00183 
00184     ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
00185     IFPACK_CHK_ERR(ierr);
00186 
00187     for (int j = 0 ; j < NumIndices ; ++j) {
00188       int jj = Indices[j];
00189       if (jj != i) {
00190     adjncy[count++] = jj;
00191     xadj[count2+1]++;
00192       }
00193     }
00194     count2++;
00195   }
00196 
00197   std::vector<idxtype> NodesInSubgraph;
00198   NodesInSubgraph.resize(NumLocalParts_);
00199 
00200   // some cases can be handled separately
00201   
00202   int ok;
00203 
00204   if (NumLocalParts() == 1) {
00205 
00206     for (int i = 0 ; i < NumMyRows() ; ++i) 
00207       Partition_[i] = 0;
00208     
00209   } else if (NumLocalParts() == NumMyRows()) {
00210 
00211     for (int i = 0 ; i < NumMyRows() ; ++i) 
00212       Partition_[i] = i;
00213   
00214   } else {
00215 
00216     ok = 0;
00217 
00218     // sometimes METIS creates less partitions than specified.
00219     // ok will check this problem, and recall metis, asking
00220     // for NumLocalParts_/2 partitions
00221     while (ok == 0) {
00222       
00223       for (int i = 0 ; i < NumMyRows() ; ++i) 
00224     Partition_[i] = -1;
00225     
00226 #ifdef HAVE_IFPACK_METIS
00227       int j = NumMyRows();
00228       if (NumLocalParts_ < 8) {
00229 
00230     int i = 1; /* optype in the METIS manual */
00231     numflag = 0;
00232     METIS_EstimateMemory(&j, &xadj[0], &adjncy[0], 
00233                  &numflag, &i, &nbytes );
00234     
00235     METIS_PartGraphRecursive(&j, &xadj[0], &adjncy[0],
00236                  NULL, NULL,
00237                  &wgtflag[0], &numflag, &NumLocalParts_, 
00238                  &options[0], &edgecut, &Partition_[0]);
00239       } else {
00240 
00241     numflag = 0;
00242     
00243     METIS_PartGraphKway (&j, &xadj[0], &adjncy[0], 
00244                  NULL, 
00245                  NULL, &wgtflag[0], &numflag, 
00246                  &NumLocalParts_, &options[0],
00247                  &edgecut, &Partition_[0]);
00248       }
00249 #else
00250       numflag = numflag * 2; // avoid warning for unused variable
00251       if (Graph_->Comm().MyPID() == 0) {
00252     cerr << "METIS was not linked; now I put all" << endl;
00253     cerr << "the local nodes in the same partition." << endl;
00254       }
00255       for (int i = 0 ; i < NumMyRows() ; ++i) 
00256     Partition_[i] = 0;
00257       NumLocalParts_ = 1;
00258 #endif
00259       
00260       ok = 1;
00261       
00262       for (int i = 0 ; i < NumLocalParts() ; ++i) 
00263     NodesInSubgraph[i] = 0;
00264 
00265       for (int i = 0 ; i < NumMyRows() ; ++i) {
00266     int j = Partition_[i];
00267     if ((j < 0) || (j>= NumLocalParts())) {
00268       ok = 0;
00269       break;
00270     } 
00271     else NodesInSubgraph[j]++;
00272       }
00273       
00274       for (int i = 0 ; i < NumLocalParts() ; ++i) {
00275     if( NodesInSubgraph[i] == 0 ) {
00276       ok = 0;
00277       break;
00278     }
00279       }
00280       
00281       if (ok == 0) {
00282     cerr << "Specified number of subgraphs ("
00283          << NumLocalParts_ << ") generates empty subgraphs." << endl;
00284     cerr << "Now I recall METIS with NumLocalParts_ = "
00285          << NumLocalParts_ / 2 << "..." << endl;
00286     NumLocalParts_ = NumLocalParts_/2;
00287       }
00288       
00289       if (NumLocalParts() == 0) {
00290     IFPACK_CHK_ERR(-10); // something went wrong
00291       }
00292       
00293       if (NumLocalParts() == 1) {
00294     for (int i = 0 ; i < NumMyRows() ; ++i) 
00295       Partition_[i] = 0;
00296     ok = 1;
00297       }
00298       
00299     } /* while( ok == 0 ) */
00300   
00301   } /* if( NumLocalParts_ == 1 ) */
00302 
00303   return(0);
00304 } 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends