Ifpack Package Browser (Single Doxygen Collection)  Development
Ifpack_AMDReordering.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 "Teuchos_ParameterList.hpp"
45 #include "Teuchos_RefCountPtr.hpp"
46 #include "Epetra_MultiVector.h"
47 #include "Ifpack_Graph.h"
48 #include "Epetra_RowMatrix.h"
50 #include "Ifpack_AMDReordering.h"
51 
52 extern "C" {
53 #include <amesos_amd.h>
54 }
55 
56 //==============================================================================
59  NumMyRows_(0),
60  IsComputed_(false)
61 {
62 }
63 
64 //==============================================================================
67  NumMyRows_(RHS.NumMyRows()),
68  IsComputed_(RHS.IsComputed())
69 {
70  Reorder_.resize(NumMyRows());
71  InvReorder_.resize(NumMyRows());
72  for (int i = 0 ; i < NumMyRows() ; ++i) {
73  Reorder_[i] = RHS.Reorder(i);
74  InvReorder_[i] = RHS.InvReorder(i);
75  }
76 }
77 
78 //==============================================================================
81 {
82  if (this == &RHS) {
83  return (*this);
84  }
85 
86  NumMyRows_ = RHS.NumMyRows(); // set number of local rows
87  IsComputed_ = RHS.IsComputed();
88  // resize vectors, and copy values from RHS
89  Reorder_.resize(NumMyRows());
90  InvReorder_.resize(NumMyRows());
91  if (IsComputed()) {
92  for (int i = 0 ; i < NumMyRows_ ; ++i) {
93  Reorder_[i] = RHS.Reorder(i);
94  InvReorder_[i] = RHS.InvReorder(i);
95  }
96  }
97  return (*this);
98 }
99 
100 //==============================================================================
102 SetParameter(const std::string Name, const int Value)
103 {
104  return(0);
105 }
106 
107 //==============================================================================
109 SetParameter(const std::string Name, const double Value)
110 {
111  return(0);
112 }
113 
114 //==============================================================================
117 {
118  return(0);
119 }
120 
121 //==============================================================================
123 {
125 
127 
128  return(0);
129 }
130 
131 //==============================================================================
133 {
134  using std::cout;
135  using std::endl;
136 
137  IsComputed_ = false;
138  NumMyRows_ = Graph.NumMyRows();
139  int NumNz = Graph.NumMyNonzeros();
140 
141  if (NumMyRows_ == 0)
142  IFPACK_CHK_ERR(-1); // strange graph this one
143 
144  // Extract CRS format
145  std::vector<int> ia(NumMyRows_+1,0);
146  std::vector<int> ja(NumNz);
147  int cnt;
148  for( int i = 0; i < NumMyRows_; ++i )
149  {
150  int * tmpP = &ja[ia[i]];
151  Graph.ExtractMyRowCopy( i, NumNz-ia[i], cnt, tmpP );
152  ia[i+1] = ia[i] + cnt;
153  }
154 
155  // Trim down to local only
156  std::vector<int> iat(NumMyRows_+1);
157  std::vector<int> jat(NumNz);
158  int loc = 0;
159  for( int i = 0; i < NumMyRows_; ++i )
160  {
161  iat[i] = loc;
162  for( int j = ia[i]; j < ia[i+1]; ++j )
163  {
164  if( ja[j] < NumMyRows_ )
165  jat[loc++] = ja[j];
166  else
167  break;
168  }
169  }
170  iat[NumMyRows_] = loc;
171 
172  // Compute AMD permutation
173  Reorder_.resize(NumMyRows_);
174  std::vector<double> info(AMD_INFO);
175 
176  amesos_amd_order( NumMyRows_, &iat[0], &jat[0], &Reorder_[0], NULL, &info[0] );
177 
178  if( info[AMD_STATUS] == AMD_INVALID )
179  cout << "AMD ORDERING: Invalid!!!!" << endl;
180 
181  // Build inverse reorder (will be used by ExtractMyRowCopy()
182  InvReorder_.resize(NumMyRows_);
183 
184  for (int i = 0 ; i < NumMyRows_ ; ++i)
185  InvReorder_[i] = -1;
186 
187  for (int i = 0 ; i < NumMyRows_ ; ++i)
188  InvReorder_[Reorder_[i]] = i;
189 
190  for (int i = 0 ; i < NumMyRows_ ; ++i) {
191  if (InvReorder_[i] == -1)
192  IFPACK_CHK_ERR(-1);
193  }
194 
195  IsComputed_ = true;
196  return(0);
197 }
198 
199 //==============================================================================
200 int Ifpack_AMDReordering::Reorder(const int i) const
201 {
202 #ifdef IFPACK_ABC
203  if (!IsComputed())
204  IFPACK_CHK_ERR(-1);
205  if ((i < 0) || (i >= NumMyRows_))
206  IFPACK_CHK_ERR(-1);
207 #endif
208 
209  return(Reorder_[i]);
210 }
211 
212 //==============================================================================
213 int Ifpack_AMDReordering::InvReorder(const int i) const
214 {
215 #ifdef IFPACK_ABC
216  if (!IsComputed())
217  IFPACK_CHK_ERR(-1);
218  if ((i < 0) || (i >= NumMyRows_))
219  IFPACK_CHK_ERR(-1);
220 #endif
221 
222  return(InvReorder_[i]);
223 }
224 //==============================================================================
226  Epetra_MultiVector& X) const
227 {
228  int NumVectors = X.NumVectors();
229 
230  for (int j = 0 ; j < NumVectors ; ++j) {
231  for (int i = 0 ; i < NumMyRows_ ; ++i) {
232  int np = Reorder_[i];
233  X[j][np] = Xorig[j][i];
234  }
235  }
236 
237  return(0);
238 }
239 
240 //==============================================================================
242  Epetra_MultiVector& X) const
243 {
244  int NumVectors = X.NumVectors();
245 
246  for (int j = 0 ; j < NumVectors ; ++j) {
247  for (int i = 0 ; i < NumMyRows_ ; ++i) {
248  int np = Reorder_[i];
249  X[j][i] = Xorig[j][np];
250  }
251  }
252 
253  return(0);
254 }
255 
256 //==============================================================================
257 std::ostream& Ifpack_AMDReordering::Print(std::ostream& os) const
258 {
259  using std::endl;
260 
261  os << "*** Ifpack_AMDReordering" << endl << endl;
262  if (!IsComputed())
263  os << "*** Reordering not yet computed." << endl;
264 
265  os << "*** Number of local rows = " << NumMyRows_ << endl;
266  os << endl;
267  os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
268  for (int i = 0 ; i < NumMyRows_ ; ++i) {
269  os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
270  }
271 
272  return(os);
273 }
int NumMyRows_
Number of local rows in the graph.
std::vector< int > Reorder_
Contains the reordering.
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
int InvReorder(const int i) const
Returns the inverse reordered index of row i.
Ifpack_AMDReordering()
Constructor for Ifpack_Graph&#39;s.
int SetParameters(Teuchos::ParameterList &List)
Sets all parameters.
Ifpack_AMDReordering: approximate minimum degree reordering.
int Pinv(const Epetra_MultiVector &Xorig, Epetra_MultiVector &Xinvreord) const
Applies inverse reordering to multivector X, whose local length equals the number of local rows...
Ifpack_AMDReordering & operator=(const Ifpack_AMDReordering &RHS)
Assignment operator.
const int NumVectors
Definition: performance.cpp:71
#define AMD_STATUS
int NumMyRows() const
Returns the number of local rows.
int Compute(const Ifpack_Graph &Graph)
Computes all it is necessary to initialize the reordering object.
int SetParameter(const std::string Name, const int Value)
Sets integer parameters ‘Name’.
#define NULL
bool IsComputed() const
Returns true is the reordering object has been successfully initialized, false otherwise.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
std::vector< int > InvReorder_
Contains the inverse reordering.
int NumVectors() const
#define AMD_INVALID
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
int amesos_amd_order(int n, const int Ap [], const int Ai [], int P [], double Control [], double Info [])
int Reorder(const int i) const
Returns the reordered index of row i.
Ifpack_Graph: a pure virtual class that defines graphs for IFPACK.
Definition: Ifpack_Graph.h:61
int P(const Epetra_MultiVector &Xorig, Epetra_MultiVector &Xreord) const
Applies reordering to multivector X, whose local length equals the number of local rows...
#define IFPACK_CHK_ERR(ifpack_err)
bool IsComputed_
If true, the reordering has been successfully computed.
#define AMD_INFO
#define RHS(a)
Definition: MatGenFD.c:60