IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_METISPartitioner.cpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_Partitioner.h"
45 #include "Ifpack_OverlappingPartitioner.h"
46 #include "Ifpack_METISPartitioner.h"
47 #include "Ifpack_Graph.h"
48 #include "Ifpack_Graph_Epetra_CrsGraph.h"
49 #include "Epetra_Comm.h"
50 #include "Epetra_CrsGraph.h"
51 #include "Epetra_Map.h"
52 #include "Teuchos_ParameterList.hpp"
53 #include "Teuchos_RefCountPtr.hpp"
54 
55 // may need to change this for wierd installations
56 typedef int idxtype;
57 #ifdef HAVE_IFPACK_METIS
58 extern "C" {
59  void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *,
60  idxtype *, int *, int *, int *, int *, int *,
61  idxtype *);
62  void METIS_PartGraphRecursive(int *, idxtype *, idxtype *,
63  idxtype *, idxtype *, int *, int *, int *,
64  int *, int *, idxtype *);
65 
66 }
67 #endif
68 
69 //==============================================================================
70 // NOTE:
71 // - matrix is supposed to be localized, and passes through the
72 // singleton filter. This means that I do not have to look
73 // for Dirichlet nodes (singletons). Also, all rows and columns are
74 // local.
76 {
77  using std::cerr;
78  using std::endl;
79 
80  int ierr;
81 #ifdef HAVE_IFPACK_METIS
82  int edgecut;
83 #endif
84 
85  Teuchos::RefCountPtr<Epetra_CrsGraph> SymGraph ;
86  Teuchos::RefCountPtr<Epetra_Map> SymMap;
87  Teuchos::RefCountPtr<Ifpack_Graph_Epetra_CrsGraph> SymIFPACKGraph;
88  Teuchos::RefCountPtr<Ifpack_Graph> IFPACKGraph = Teuchos::rcp( (Ifpack_Graph*)Graph_, false );
89 
90  int Length = 2 * MaxNumEntries();
91  int NumIndices;
92  std::vector<int> Indices;
93  Indices.resize(Length);
94 
95  /* construct the CSR graph information of the LOCAL matrix
96  using the get_row function */
97 
98  std::vector<idxtype> wgtflag;
99  wgtflag.resize(4);
100 
101  std::vector<int> options;
102  options.resize(4);
103 
104  int numflag;
105 
106  if (UseSymmetricGraph_) {
107 
108 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
109  // need to build a symmetric graph.
110  // I do this in two stages:
111  // 1.- construct an Epetra_CrsMatrix, symmetric
112  // 2.- convert the Epetra_CrsMatrix into METIS format
113  SymMap = Teuchos::rcp( new Epetra_Map(NumMyRows(),0,Graph_->Comm()) );
114  SymGraph = Teuchos::rcp( new Epetra_CrsGraph(Copy,*SymMap,0) );
115 #endif
116 
117 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
118  if(SymGraph->RowMap().GlobalIndicesInt()) {
119  for (int i = 0; i < NumMyRows() ; ++i) {
120 
121  ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
122  IFPACK_CHK_ERR(ierr);
123 
124  for (int j = 0 ; j < NumIndices ; ++j) {
125  int jj = Indices[j];
126  if (jj != i) {
127  SymGraph->InsertGlobalIndices(i,1,&jj);
128  SymGraph->InsertGlobalIndices(jj,1,&i);
129  }
130  }
131  }
132  }
133  else
134 #endif
135 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
136  if(SymGraph->RowMap().GlobalIndicesLongLong()) {
137  for (int i = 0; i < NumMyRows() ; ++i) {
138  long long i_LL = i;
139 
140  ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
141  IFPACK_CHK_ERR(ierr);
142 
143  for (int j = 0 ; j < NumIndices ; ++j) {
144  long long jj = Indices[j];
145  if (jj != i_LL) {
146  SymGraph->InsertGlobalIndices(i_LL,1,&jj);
147  SymGraph->InsertGlobalIndices(jj,1,&i_LL);
148  }
149  }
150  }
151  }
152  else
153 #endif
154  throw "Ifpack_METISPartitioner::ComputePartitions: GlobalIndices type unknown";
155 
156  IFPACK_CHK_ERR(SymGraph->FillComplete());
157  SymIFPACKGraph = Teuchos::rcp( new Ifpack_Graph_Epetra_CrsGraph(SymGraph) );
158  IFPACKGraph = SymIFPACKGraph;
159  }
160 
161  // now work on IFPACKGraph, that can be the symmetric or
162  // the non-symmetric one
163 
164  /* set parameters */
165 
166  wgtflag[0] = 0; /* no weights */
167  numflag = 0; /* C style */
168  options[0] = 0; /* default options */
169 
170  std::vector<idxtype> xadj;
171  xadj.resize(NumMyRows() + 1);
172 
173  std::vector<idxtype> adjncy;
174  adjncy.resize(NumMyNonzeros());
175 
176  int count = 0;
177  int count2 = 0;
178  xadj[0] = 0;
179 
180  for (int i = 0; i < NumMyRows() ; ++i) {
181 
182  xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
183 
184  ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
185  IFPACK_CHK_ERR(ierr);
186 
187  for (int j = 0 ; j < NumIndices ; ++j) {
188  int jj = Indices[j];
189  if (jj != i) {
190  adjncy[count++] = jj;
191  xadj[count2+1]++;
192  }
193  }
194  count2++;
195  }
196 
197  std::vector<idxtype> NodesInSubgraph;
198  NodesInSubgraph.resize(NumLocalParts_);
199 
200  // some cases can be handled separately
201 
202  int ok;
203 
204  if (NumLocalParts() == 1) {
205 
206  for (int i = 0 ; i < NumMyRows() ; ++i)
207  Partition_[i] = 0;
208 
209  } else if (NumLocalParts() == NumMyRows()) {
210 
211  for (int i = 0 ; i < NumMyRows() ; ++i)
212  Partition_[i] = i;
213 
214  } else {
215 
216  ok = 0;
217 
218  // sometimes METIS creates less partitions than specified.
219  // ok will check this problem, and recall metis, asking
220  // for NumLocalParts_/2 partitions
221  while (ok == 0) {
222 
223  for (int i = 0 ; i < NumMyRows() ; ++i)
224  Partition_[i] = -1;
225 
226 #ifdef HAVE_IFPACK_METIS
227  int j = NumMyRows();
228  if (NumLocalParts_ < 8) {
229 
230  numflag = 0;
231 
232  METIS_PartGraphRecursive(&j, &xadj[0], &adjncy[0],
233  NULL, NULL,
234  &wgtflag[0], &numflag, &NumLocalParts_,
235  &options[0], &edgecut, &Partition_[0]);
236  } else {
237 
238  numflag = 0;
239 
240  METIS_PartGraphKway (&j, &xadj[0], &adjncy[0],
241  NULL,
242  NULL, &wgtflag[0], &numflag,
243  &NumLocalParts_, &options[0],
244  &edgecut, &Partition_[0]);
245  }
246 #else
247  numflag = numflag * 2; // avoid warning for unused variable
248  if (Graph_->Comm().MyPID() == 0) {
249  cerr << "METIS was not linked; now I put all" << endl;
250  cerr << "the local nodes in the same partition." << endl;
251  }
252  for (int i = 0 ; i < NumMyRows() ; ++i)
253  Partition_[i] = 0;
254  NumLocalParts_ = 1;
255 #endif
256 
257  ok = 1;
258 
259  for (int i = 0 ; i < NumLocalParts() ; ++i)
260  NodesInSubgraph[i] = 0;
261 
262  for (int i = 0 ; i < NumMyRows() ; ++i) {
263  int j = Partition_[i];
264  if ((j < 0) || (j>= NumLocalParts())) {
265  ok = 0;
266  break;
267  }
268  else NodesInSubgraph[j]++;
269  }
270 
271  for (int i = 0 ; i < NumLocalParts() ; ++i) {
272  if( NodesInSubgraph[i] == 0 ) {
273  ok = 0;
274  break;
275  }
276  }
277 
278  if (ok == 0) {
279  cerr << "Specified number of subgraphs ("
280  << NumLocalParts_ << ") generates empty subgraphs." << endl;
281  cerr << "Now I recall METIS with NumLocalParts_ = "
282  << NumLocalParts_ / 2 << "..." << endl;
284  }
285 
286  if (NumLocalParts() == 0) {
287  IFPACK_CHK_ERR(-10); // something went wrong
288  }
289 
290  if (NumLocalParts() == 1) {
291  for (int i = 0 ; i < NumMyRows() ; ++i)
292  Partition_[i] = 0;
293  ok = 1;
294  }
295 
296  } /* while( ok == 0 ) */
297 
298  } /* if( NumLocalParts_ == 1 ) */
299 
300  return(0);
301 }
virtual int ExtractMyRowCopy(int MyRow, int LenOfIndices, int &NumIndices, int *Indices) const =0
Extracts a copy of input local row.
int MaxNumEntries() const
Returns the max number of local entries in a row.
int NumMyNonzeros() const
Returns the number of local nonzero elements.
int NumMyRows() const
Returns the number of local rows.
std::vector< int > Partition_
Partition_[i] contains the ID of non-overlapping part it belongs to.
virtual int MyPID() const =0
Ifpack_Graph_Epetra_CrsGraph: a class to define Ifpack_Graph as a light-weight conversion of Epetra_C...
virtual const Epetra_Comm & Comm() const =0
Returns the communicator object of the graph.
const Ifpack_Graph * Graph_
Reference to the graph to be partitioned.
int NumLocalParts_
Number of local subgraphs.
int NumLocalParts() const
Returns the number of computed local partitions.
Ifpack_Graph: a pure virtual class that defines graphs for IFPACK.
Definition: Ifpack_Graph.h:61
int ComputePartitions()
Computes the partitions. Returns 0 if successful.