IFPACK Development
Loading...
Searching...
No Matches
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
56typedef int idxtype;
57#ifdef HAVE_IFPACK_METIS
58extern "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}
Copy
virtual int MyPID() const=0
Ifpack_Graph_Epetra_CrsGraph: a class to define Ifpack_Graph as a light-weight conversion of Epetra_C...
Ifpack_Graph: a pure virtual class that defines graphs for IFPACK.
virtual int ExtractMyRowCopy(int MyRow, int LenOfIndices, int &NumIndices, int *Indices) const =0
Extracts a copy of input local row.
virtual const Epetra_Comm & Comm() const =0
Returns the communicator object of the graph.
int ComputePartitions()
Computes the partitions. Returns 0 if successful.
int NumLocalParts() const
Returns the number of computed local partitions.
int NumMyNonzeros() const
Returns the number of local nonzero elements.
int NumLocalParts_
Number of local subgraphs.
int NumMyRows() const
Returns the number of local rows.
int MaxNumEntries() const
Returns the max number of local entries in a row.
std::vector< int > Partition_
Partition_[i] contains the ID of non-overlapping part it belongs to.
const Ifpack_Graph * Graph_
Reference to the graph to be partitioned.