IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_AMDReordering.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 "Teuchos_ParameterList.hpp"
45 #include "Teuchos_RefCountPtr.hpp"
46 #include "Epetra_MultiVector.h"
47 #include "Ifpack_Graph.h"
48 #include "Epetra_RowMatrix.h"
49 #include "Ifpack_Graph_Epetra_RowMatrix.h"
50 #include "Ifpack_AMDReordering.h"
51 
52 extern "C" {
53 #include <trilinos_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 //==============================================================================
116 SetParameters(Teuchos::ParameterList& /* List */)
117 {
118  return(0);
119 }
120 
121 //==============================================================================
123 {
124  Ifpack_Graph_Epetra_RowMatrix Graph(Teuchos::rcp(&Matrix,false));
125 
126  IFPACK_CHK_ERR(Compute(Graph));
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(TRILINOS_AMD_INFO);
175 
176  trilinos_amd_order( NumMyRows_, &iat[0], &jat[0], &Reorder_[0], NULL, &info[0] );
177 
178  if( info[TRILINOS_AMD_STATUS] == TRILINOS_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 }
virtual int ExtractMyRowCopy(int MyRow, int LenOfIndices, int &NumIndices, int *Indices) const =0
Extracts a copy of input local row.
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...
virtual int NumMyRows() const =0
Returns the number of local rows.
Ifpack_AMDReordering()
Constructor for Ifpack_Graph&#39;s.
int SetParameters(Teuchos::ParameterList &List)
Sets all parameters.
Ifpack_AMDReordering: approximate minimum degree reordering.
Ifpack_AMDReordering & operator=(const Ifpack_AMDReordering &RHS)
Assignment operator.
int P(const Epetra_MultiVector &Xorig, Epetra_MultiVector &Xreord) const
Applies reordering to multivector X, whose local length equals the number of local rows...
bool IsComputed() const
Returns true is the reordering object has been successfully initialized, false otherwise.
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&#39;.
virtual int NumMyNonzeros() const =0
Returns the number of local nonzero entries.
int Reorder(const int i) const
Returns the reordered index of row i.
int InvReorder(const int i) const
Returns the inverse reordered index of row i.
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
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator&lt;&lt;.