Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
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
55typedef int idxtype;
56#ifdef HAVE_IFPACK_METIS
57extern "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//==============================================================================
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//==============================================================================
283std::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
Copy
#define IFPACK_CHK_ERR(ifpack_err)
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
int NumVectors() const
Ifpack_Graph_Epetra_CrsGraph: a class to define Ifpack_Graph as a light-weight conversion of Epetra_C...
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
Ifpack_Graph: a pure virtual class that defines graphs for IFPACK.
Definition: Ifpack_Graph.h:61
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.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
bool UseSymmetricGraph_
If true, the graph has to be symmetrized before calling METIS.
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,...
virtual bool IsComputed() const
Returns true is the reordering object has been successfully initialized, false otherwise.
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 int Reorder(const int i) const
Returns the reordered index of row i.
int NumMyRows_
Number of local rows in the graph.
std::vector< int > InvReorder_
Contains the inverse reordering.
#define false
const int NumVectors
Definition: performance.cpp:71