Ifpack Package Browser (Single Doxygen Collection)  Development
Ifpack_METISPartitioner.cpp
Go to the documentation of this file.
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"
47 #include "Ifpack_Graph.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_EstimateMemory(int *, idxtype *, idxtype *, int *, int *, int *);
60  void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *,
61  idxtype *, int *, int *, int *, int *, int *,
62  idxtype *);
63  void METIS_PartGraphRecursive(int *, idxtype *, idxtype *,
64  idxtype *, idxtype *, int *, int *, int *,
65  int *, int *, idxtype *);
66 
67 }
68 #endif
69 
70 //==============================================================================
71 // NOTE:
72 // - matrix is supposed to be localized, and passes through the
73 // singleton filter. This means that I do not have to look
74 // for Dirichlet nodes (singletons). Also, all rows and columns are
75 // local.
77 {
78  using std::cerr;
79  using std::endl;
80 
81  int ierr;
82 #ifdef HAVE_IFPACK_METIS
83  int nbytes = 0;
84  int edgecut;
85 #endif
86 
87  Teuchos::RefCountPtr<Epetra_CrsGraph> SymGraph ;
88  Teuchos::RefCountPtr<Epetra_Map> SymMap;
89  Teuchos::RefCountPtr<Ifpack_Graph_Epetra_CrsGraph> SymIFPACKGraph;
90  Teuchos::RefCountPtr<Ifpack_Graph> IFPACKGraph = Teuchos::rcp( (Ifpack_Graph*)Graph_, false );
91 
92  int Length = 2 * MaxNumEntries();
93  int NumIndices;
94  std::vector<int> Indices;
95  Indices.resize(Length);
96 
97  /* construct the CSR graph information of the LOCAL matrix
98  using the get_row function */
99 
100  std::vector<idxtype> wgtflag;
101  wgtflag.resize(4);
102 
103  std::vector<int> options;
104  options.resize(4);
105 
106  int numflag;
107 
108  if (UseSymmetricGraph_) {
109 
110 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
111  // need to build a symmetric graph.
112  // I do this in two stages:
113  // 1.- construct an Epetra_CrsMatrix, symmetric
114  // 2.- convert the Epetra_CrsMatrix into METIS format
115  SymMap = Teuchos::rcp( new Epetra_Map(NumMyRows(),0,Graph_->Comm()) );
116  SymGraph = Teuchos::rcp( new Epetra_CrsGraph(Copy,*SymMap,0) );
117 #endif
118 
119 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
120  if(SymGraph->RowMap().GlobalIndicesInt()) {
121  for (int i = 0; i < NumMyRows() ; ++i) {
122 
123  ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
124  IFPACK_CHK_ERR(ierr);
125 
126  for (int j = 0 ; j < NumIndices ; ++j) {
127  int jj = Indices[j];
128  if (jj != i) {
129  SymGraph->InsertGlobalIndices(i,1,&jj);
130  SymGraph->InsertGlobalIndices(jj,1,&i);
131  }
132  }
133  }
134  }
135  else
136 #endif
137 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
138  if(SymGraph->RowMap().GlobalIndicesLongLong()) {
139  for (int i = 0; i < NumMyRows() ; ++i) {
140  long long i_LL = i;
141 
142  ierr = Graph_->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
143  IFPACK_CHK_ERR(ierr);
144 
145  for (int j = 0 ; j < NumIndices ; ++j) {
146  long long jj = Indices[j];
147  if (jj != i_LL) {
148  SymGraph->InsertGlobalIndices(i_LL,1,&jj);
149  SymGraph->InsertGlobalIndices(jj,1,&i_LL);
150  }
151  }
152  }
153  }
154  else
155 #endif
156  throw "Ifpack_METISPartitioner::ComputePartitions: GlobalIndices type unknown";
157 
158  IFPACK_CHK_ERR(SymGraph->FillComplete());
159  SymIFPACKGraph = Teuchos::rcp( new Ifpack_Graph_Epetra_CrsGraph(SymGraph) );
160  IFPACKGraph = SymIFPACKGraph;
161  }
162 
163  // now work on IFPACKGraph, that can be the symmetric or
164  // the non-symmetric one
165 
166  /* set parameters */
167 
168  wgtflag[0] = 0; /* no weights */
169  numflag = 0; /* C style */
170  options[0] = 0; /* default options */
171 
172  std::vector<idxtype> xadj;
173  xadj.resize(NumMyRows() + 1);
174 
175  std::vector<idxtype> adjncy;
176  adjncy.resize(NumMyNonzeros());
177 
178  int count = 0;
179  int count2 = 0;
180  xadj[0] = 0;
181 
182  for (int i = 0; i < NumMyRows() ; ++i) {
183 
184  xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
185 
186  ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
187  IFPACK_CHK_ERR(ierr);
188 
189  for (int j = 0 ; j < NumIndices ; ++j) {
190  int jj = Indices[j];
191  if (jj != i) {
192  adjncy[count++] = jj;
193  xadj[count2+1]++;
194  }
195  }
196  count2++;
197  }
198 
199  std::vector<idxtype> NodesInSubgraph;
200  NodesInSubgraph.resize(NumLocalParts_);
201 
202  // some cases can be handled separately
203 
204  int ok;
205 
206  if (NumLocalParts() == 1) {
207 
208  for (int i = 0 ; i < NumMyRows() ; ++i)
209  Partition_[i] = 0;
210 
211  } else if (NumLocalParts() == NumMyRows()) {
212 
213  for (int i = 0 ; i < NumMyRows() ; ++i)
214  Partition_[i] = i;
215 
216  } else {
217 
218  ok = 0;
219 
220  // sometimes METIS creates less partitions than specified.
221  // ok will check this problem, and recall metis, asking
222  // for NumLocalParts_/2 partitions
223  while (ok == 0) {
224 
225  for (int i = 0 ; i < NumMyRows() ; ++i)
226  Partition_[i] = -1;
227 
228 #ifdef HAVE_IFPACK_METIS
229  int j = NumMyRows();
230  if (NumLocalParts_ < 8) {
231 
232  int i = 1; /* optype in the METIS manual */
233  numflag = 0;
234  METIS_EstimateMemory(&j, &xadj[0], &adjncy[0],
235  &numflag, &i, &nbytes );
236 
237  METIS_PartGraphRecursive(&j, &xadj[0], &adjncy[0],
238  NULL, NULL,
239  &wgtflag[0], &numflag, &NumLocalParts_,
240  &options[0], &edgecut, &Partition_[0]);
241  } else {
242 
243  numflag = 0;
244 
245  METIS_PartGraphKway (&j, &xadj[0], &adjncy[0],
246  NULL,
247  NULL, &wgtflag[0], &numflag,
248  &NumLocalParts_, &options[0],
249  &edgecut, &Partition_[0]);
250  }
251 #else
252  numflag = numflag * 2; // avoid warning for unused variable
253  if (Graph_->Comm().MyPID() == 0) {
254  cerr << "METIS was not linked; now I put all" << endl;
255  cerr << "the local nodes in the same partition." << endl;
256  }
257  for (int i = 0 ; i < NumMyRows() ; ++i)
258  Partition_[i] = 0;
259  NumLocalParts_ = 1;
260 #endif
261 
262  ok = 1;
263 
264  for (int i = 0 ; i < NumLocalParts() ; ++i)
265  NodesInSubgraph[i] = 0;
266 
267  for (int i = 0 ; i < NumMyRows() ; ++i) {
268  int j = Partition_[i];
269  if ((j < 0) || (j>= NumLocalParts())) {
270  ok = 0;
271  break;
272  }
273  else NodesInSubgraph[j]++;
274  }
275 
276  for (int i = 0 ; i < NumLocalParts() ; ++i) {
277  if( NodesInSubgraph[i] == 0 ) {
278  ok = 0;
279  break;
280  }
281  }
282 
283  if (ok == 0) {
284  cerr << "Specified number of subgraphs ("
285  << NumLocalParts_ << ") generates empty subgraphs." << endl;
286  cerr << "Now I recall METIS with NumLocalParts_ = "
287  << NumLocalParts_ / 2 << "..." << endl;
289  }
290 
291  if (NumLocalParts() == 0) {
292  IFPACK_CHK_ERR(-10); // something went wrong
293  }
294 
295  if (NumLocalParts() == 1) {
296  for (int i = 0 ; i < NumMyRows() ; ++i)
297  Partition_[i] = 0;
298  ok = 1;
299  }
300 
301  } /* while( ok == 0 ) */
302 
303  } /* if( NumLocalParts_ == 1 ) */
304 
305  return(0);
306 }
virtual int ExtractMyRowCopy(int MyRow, int LenOfIndices, int &NumIndices, int *Indices) const =0
Extracts a copy of input local row.
Copy
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.
int NumMyRows() const
Returns the number of local rows.
int NumLocalParts() const
Returns the number of computed local partitions.
const Ifpack_Graph * Graph_
Reference to the graph to be partitioned.
int NumLocalParts_
Number of local subgraphs.
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.
#define IFPACK_CHK_ERR(ifpack_err)
int MaxNumEntries() const
Returns the max number of local entries in a row.
int NumMyNonzeros() const
Returns the number of local nonzero elements.