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