Ifpack Package Browser (Single Doxygen Collection)  Development
Ifpack_METISReordering.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_Reordering.h"
45 #include "Ifpack_METISReordering.h"
46 #include "Ifpack_Graph.h"
49 #include "Epetra_Comm.h"
50 #include "Epetra_MultiVector.h"
51 #include "Epetra_CrsGraph.h"
52 #include "Epetra_Map.h"
53 #include "Teuchos_ParameterList.hpp"
54 
55 typedef int idxtype;
56 #ifdef HAVE_IFPACK_METIS
57 extern "C" {
58  void METIS_NodeND(int *n, idxtype *xadj, idxtype *adjncy,
59  int *numflag, int *options, int *perm, int *iperm);
60 }
61 #endif
62 
63 //==============================================================================
65  UseSymmetricGraph_(false),
66  NumMyRows_(0),
67  IsComputed_(false)
68 {}
69 
70 //==============================================================================
71 // Mainly copied from Ifpack_METISPartitioner.cpp
72 //
73 // NOTE:
74 // - matrix is supposed to be localized, and passes through the
75 // singleton filter. This means that I do not have to look
76 // for Dirichlet nodes (singletons). Also, all rows and columns are
77 // local.
79 {
80 
81  NumMyRows_ = Graph.NumMyRows();
82  Reorder_.resize(NumMyRows_);
83  InvReorder_.resize(NumMyRows_);
84 
85  int ierr;
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 * Graph.MaxMyNumEntries();
93  int NumIndices;
94  std::vector<int> Indices;
95  Indices.resize(Length);
96 
97  std::vector<int> options;
98  options.resize(8);
99  options[0] = 0; // default values
100 
101 #ifdef HAVE_IFPACK_METIS
102  int numflag = 0; // C style
103 #endif
104 
105  if (UseSymmetricGraph_) {
106 
107 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
108  // need to build a symmetric graph.
109  // I do this in two stages:
110  // 1.- construct an Epetra_CrsMatrix, symmetric
111  // 2.- convert the Epetra_CrsMatrix into METIS format
112  SymMap = Teuchos::rcp( new Epetra_Map(NumMyRows_,0,Graph.Comm()) );
113  SymGraph = Teuchos::rcp( new Epetra_CrsGraph(Copy,*SymMap,0) );
114 #endif
115 
116 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
117  if(SymGraph->RowMap().GlobalIndicesInt()) {
118  for (int i = 0; i < NumMyRows_ ; ++i) {
119 
120  ierr = Graph.ExtractMyRowCopy(i, Length, NumIndices,
121  &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  // insert A(i,j), then A(j,i)
128  SymGraph->InsertGlobalIndices(i,1,&jj);
129  SymGraph->InsertGlobalIndices(jj,1,&i);
130  }
131  }
132  }
133  }
134  else
135 #endif
136 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
137  if(SymGraph->RowMap().GlobalIndicesLongLong()) {
138  for (int i = 0; i < NumMyRows_ ; ++i) {
139  long long i_LL = i;
140 
141  ierr = Graph.ExtractMyRowCopy(i, Length, NumIndices,
142  &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) {
148  // insert A(i,j), then A(j,i)
149  SymGraph->InsertGlobalIndices(i_LL,1,&jj);
150  SymGraph->InsertGlobalIndices(jj,1,&i_LL);
151  }
152  }
153  }
154  }
155  else
156 #endif
157  throw "Ifpack_METISReordering::Compute: GlobalIndices type unknown";
158 
159  IFPACK_CHK_ERR(SymGraph->OptimizeStorage());
160  IFPACK_CHK_ERR(SymGraph->FillComplete());
161  SymIFPACKGraph = Teuchos::rcp( new Ifpack_Graph_Epetra_CrsGraph(SymGraph) );
162  IFPACKGraph = SymIFPACKGraph;
163  }
164 
165  // convert to METIS format
166  std::vector<idxtype> xadj;
167  xadj.resize(NumMyRows_ + 1);
168 
169  std::vector<idxtype> adjncy;
170  adjncy.resize(Graph.NumMyNonzeros());
171 
172  int count = 0;
173  int count2 = 0;
174  xadj[0] = 0;
175 
176  for (int i = 0; i < NumMyRows_ ; ++i) {
177 
178  xadj[count2+1] = xadj[count2]; /* nonzeros in row i-1 */
179 
180  ierr = IFPACKGraph->ExtractMyRowCopy(i, Length, NumIndices, &Indices[0]);
181  IFPACK_CHK_ERR(ierr);
182 
183  for (int j = 0 ; j < NumIndices ; ++j) {
184  int jj = Indices[j];
185  if (jj != i) {
186  adjncy[count++] = jj;
187  xadj[count2+1]++;
188  }
189  }
190  count2++;
191  }
192 
193 #ifdef HAVE_IFPACK_METIS
194  // vectors from METIS. The second last is `perm', the last is `iperm'.
195  // They store the fill-reducing permutation and inverse-permutation.
196  // Let A be the original matrix and A' the permuted matrix. The
197  // arrays perm and iperm are defined as follows. Row (column) i of A'
198  // if the perm[i] row (col) of A, and row (column) i of A is the
199  // iperm[i] row (column) of A'. The numbering starts from 0 in our case.
200  METIS_NodeND(&NumMyRows_, &xadj[0], &adjncy[0],
201  &numflag, &options[0],
202  &InvReorder_[0], &Reorder_[0]);
203 #else
204  using std::cerr;
205  using std::endl;
206 
207  cerr << "Please configure with --enable-ifpack-metis" << endl;
208  cerr << "to use METIS Reordering." << endl;
209  exit(EXIT_FAILURE);
210 #endif
211 
212  return(0);
213 }
214 
215 //==============================================================================
217 {
218  Ifpack_Graph_Epetra_RowMatrix Graph(Teuchos::rcp(&Matrix, false));
219 
221 
222  return(0);
223 }
224 
225 //==============================================================================
226 int Ifpack_METISReordering::Reorder(const int i) const
227 {
228 #ifdef IFPACK_ABC
229  if (!IsComputed())
230  IFPACK_CHK_ERR(-1);
231  if ((i < 0) || (i >= NumMyRows_))
232  IFPACK_CHK_ERR(-1);
233 #endif
234 
235  return(Reorder_[i]);
236 }
237 
238 //==============================================================================
240 {
241 #ifdef IFPACK_ABC
242  if (!IsComputed())
243  IFPACK_CHK_ERR(-1);
244  if ((i < 0) || (i >= NumMyRows_))
245  IFPACK_CHK_ERR(-1);
246 #endif
247 
248  return(InvReorder_[i]);
249 }
250 //==============================================================================
252  Epetra_MultiVector& X) const
253 {
254  int NumVectors = X.NumVectors();
255 
256  for (int j = 0 ; j < NumVectors ; ++j) {
257  for (int i = 0 ; i < NumMyRows_ ; ++i) {
258  int np = Reorder_[i];
259  X[j][np] = Xorig[j][i];
260  }
261  }
262 
263  return(0);
264 }
265 
266 //==============================================================================
268  Epetra_MultiVector& X) const
269 {
270  int NumVectors = X.NumVectors();
271 
272  for (int j = 0 ; j < NumVectors ; ++j) {
273  for (int i = 0 ; i < NumMyRows_ ; ++i) {
274  int np = Reorder_[i];
275  X[j][i] = Xorig[j][np];
276  }
277  }
278 
279  return(0);
280 }
281 
282 //==============================================================================
283 std::ostream& Ifpack_METISReordering::Print(std::ostream& os) const
284 {
285  using std::endl;
286 
287  os << "*** Ifpack_METISReordering" << endl << endl;
288  if (!IsComputed())
289  os << "*** Reordering not yet computed." << endl;
290 
291  os << "*** Number of local rows = " << NumMyRows_ << endl;
292  os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
293  for (int i = 0 ; i < NumMyRows_ ; ++i) {
294  os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
295  }
296 
297  return(os);
298 }
299 
virtual int Reorder(const int i) const
Returns the reordered index of row i.
Copy
const int NumVectors
Definition: performance.cpp:71
virtual int P(const Epetra_MultiVector &Xorig, Epetra_MultiVector &X) const
Applies reordering to multivector Xorig, whose local length equals the number of local rows...
int NumMyRows_
Number of local rows in the graph.
int NumVectors() const
Ifpack_Graph_Epetra_CrsGraph: a class to define Ifpack_Graph as a light-weight conversion of Epetra_C...
bool UseSymmetricGraph_
If true, the graph has to be symmetrized before calling METIS.
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
#define false
virtual int Compute(const Ifpack_Graph &Graph)
Computes all it is necessary to initialize the reordering object.
std::vector< int > Reorder_
Contains the reordering.
virtual bool IsComputed() const
Returns true is the reordering object has been successfully initialized, false otherwise.
Ifpack_Graph: a pure virtual class that defines graphs for IFPACK.
Definition: Ifpack_Graph.h:61
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
#define IFPACK_CHK_ERR(ifpack_err)
std::vector< int > InvReorder_
Contains the inverse reordering.
virtual int Pinv(const Epetra_MultiVector &Xorig, Epetra_MultiVector &X) const
Applies inverse reordering to multivector Xorig, whose local length equals the number of local rows...
virtual int InvReorder(const int i) const
Returns the inverse reordered index of row i.