IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_Euclid.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_Euclid.h"
44 #if defined(HAVE_EUCLID) && defined(HAVE_MPI)
45 
46 #include "Ifpack_Utils.h"
47 #include <algorithm>
48 #include "Epetra_MpiComm.h"
49 #include "Epetra_IntVector.h"
50 #include "Epetra_Import.h"
51 #include "Teuchos_ParameterList.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "getRow_dh.h"
54 
55 using Teuchos::RCP;
56 using Teuchos::rcp;
57 
58 Ifpack_Euclid::Ifpack_Euclid(Epetra_CrsMatrix* A):
59  A_(rcp(A,false)),
60  UseTranspose_(false),
61  Condest_(-1),
62  IsInitialized_(false),
63  IsComputed_(false),
64  Label_(),
65  NumInitialize_(0),
66  NumCompute_(0),
67  NumApplyInverse_(0),
68  InitializeTime_(0.0),
69  ComputeTime_(0.0),
70  ApplyInverseTime_(0.0),
71  ComputeFlops_(0.0),
72  ApplyInverseFlops_(0.0),
73  Time_(A_->Comm()),
74  SetLevel_(1),
75  SetBJ_(0),
76  SetStats_(0),
77  SetMem_(0),
78  SetSparse_(0.0),
79  SetRowScale_(0),
80  SetILUT_(0.0)
81 {
82  // Here we need to change the view of each row to have global indices. This is
83  // because Euclid directly extracts a row view and expects global indices.
84  for(int i = 0; i < A_->NumMyRows(); i++){
85  int *indices;
86  int len;
87  A_->Graph().ExtractMyRowView(i, len, indices);
88  for(int j = 0; j < len; j++){
89  indices[j] = A_->GCID(indices[j]);
90  }
91  }
92 } //Constructor
93 
94 //==============================================================================
95 void Ifpack_Euclid::Destroy(){
96  // Destroy the euclid solver, only if it was setup
97  if(IsComputed()){
98  Euclid_dhDestroy(eu);
99  }
100  // Delete these euclid varaiables if they were created
101  if(IsInitialized()){
102  Parser_dhDestroy(parser_dh);
103  parser_dh = NULL;
104  TimeLog_dhDestroy(tlog_dh);
105  tlog_dh = NULL;
106  Mem_dhDestroy(mem_dh);
107  mem_dh = NULL;
108  }
109  // Now that Euclid is done with the matrix, we change it back to having local indices.
110  for(int i = 0; i < A_->NumMyRows(); i++){
111  int *indices;
112  int len;
113  A_->Graph().ExtractMyRowView(i, len, indices);
114  for(int j = 0; j < len; j++){
115  indices[j] = A_->LCID(indices[j]);
116  }
117  }
118 } //Destroy()
119 
120 //==============================================================================
121 int Ifpack_Euclid::Initialize(){
122  //These are global variables in Euclid
123  comm_dh = GetMpiComm();
124  MPI_Comm_size(comm_dh, &np_dh);
125  MPI_Comm_rank(comm_dh, &myid_dh);
126  Time_.ResetStartTime();
127  if(mem_dh == NULL){
128  Mem_dhCreate(&mem_dh);
129  }
130  if (tlog_dh == NULL) {
131  TimeLog_dhCreate(&tlog_dh);
132  }
133 
134  if (parser_dh == NULL) {
135  Parser_dhCreate(&parser_dh);
136  }
137  Parser_dhInit(parser_dh, 0, NULL);
138  // Create the solver, this doesn't malloc anything yet, so it's only destroyed if Compute() is called.
139  Euclid_dhCreate(&eu);
140  IsInitialized_=true;
141  NumInitialize_ = NumInitialize_ + 1;
142  InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
143  return 0;
144 } //Initialize()
145 
146 //==============================================================================
147 int Ifpack_Euclid::SetParameters(Teuchos::ParameterList& list){
148  List_ = list;
149  SetLevel_ = list.get("SetLevel", (int)1);
150  SetBJ_ = list.get("SetBJ", (int)0);
151  SetStats_ = list.get("SetStats", (int)0);
152  SetMem_ = list.get("SetMem", (int)0);
153  SetSparse_ = list.get("SetSparse", (double)0.0);
154  SetRowScale_ = list.get("SetRowScale", (int)0);
155  SetILUT_ = list.get("SetILUT", (double)0.0);
156  return 0;
157 } //SetParamters()
158 
159 //==============================================================================
160 int Ifpack_Euclid::SetParameter(std::string name, int value){
161  //Convert to lowercase (so it's case insensitive)
162  locale loc;
163  for(size_t i = 0; i < name.length(); i++){
164  name[i] = (char) tolower(name[i], loc);
165  }
166  if(name.compare("setlevel") == 0){
167  SetLevel_ = value;
168  } else if(name.compare("setbj") == 0){
169  SetBJ_ = value;
170  } else if(name.compare("setstats") == 0){
171  SetStats_ = value;
172  } else if(name.compare("setmem") == 0){
173  SetMem_ = value;
174  } else if(name.compare("setrowscale") == 0){
175  SetRowScale_ = value;
176  } else {
177  using std::cout;
178  using std::endl;
179 
180  cout << "\nThe string " << name << " is not an available option." << endl;
181  IFPACK_CHK_ERR(-1);
182  }
183  return 0;
184 } //SetParameter() (int)
185 
186 //==============================================================================
187 int Ifpack_Euclid::SetParameter(std::string name, double value){
188  //Convert to lowercase (so it's case insensitive)
189  locale loc;
190  for(size_t i; i < name.length(); i++){
191  name[i] = (char) tolower(name[i], loc);
192  }
193  if(name.compare("setsparse") == 0){
194  SetSparse_ = value;
195  } else if(name.compare("setilut") == 0){
196  SetILUT_ = value;
197  } else {
198  using std::cout;
199  using std::endl;
200 
201  cout << "\nThe string " << name << " is not an available option." << endl;
202  IFPACK_CHK_ERR(-1);
203  }
204  return 0;
205 } //SetParameter() (double)
206 
207 //==============================================================================
208 int Ifpack_Euclid::Compute(){
209  if(IsInitialized() == false){
210  IFPACK_CHK_ERR(Initialize());
211  }
212  Time_.ResetStartTime();
213  sprintf(Label_, "IFPACK_Euclid (level=%d, bj=%d, stats=%d, mem=%d, sparse=%f, rowscale=%d, ilut=%f)",
214  SetLevel_, SetBJ_, SetStats_, SetMem_, SetSparse_, SetRowScale_, SetILUT_);
215  // Set the parameters
216  eu->level = SetLevel_;
217  if(SetBJ_ != 0){
218  strcpy("bj", eu->algo_par);
219  }
220  if(SetSparse_ != 0.0){
221  eu->sparseTolA = SetSparse_;
222  }
223  if(SetRowScale_ != 0){
224  eu->isScaled = true;
225  }
226  if(SetILUT_ != 0.0){
227  eu->droptol = SetILUT_;
228  }
229  if(SetStats_ != 0 || SetMem_ != 0){
230  eu->logging = true;
231  Parser_dhInsert(parser_dh, "-eu_stats", "1");
232  }
233  // eu->A is the matrix as a void pointer, eu->m is local rows, eu->n is global rows
234  eu->A = (void*) A_.get();
235  eu->m = A_->NumMyRows();
236  eu->n = A_->NumGlobalRows();
237  Euclid_dhSetup(eu);
238  IsComputed_ = true;
239  NumCompute_ = NumCompute_ + 1;
240  ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
241  return 0;
242 } //Compute()
243 
244 //==============================================================================
245 int Ifpack_Euclid::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
246  if(IsComputed() == false){
247  IFPACK_CHK_ERR(-1);
248  }
249  int NumVectors = X.NumVectors();
250  if(NumVectors != Y.NumVectors()){
251  IFPACK_CHK_ERR(-2);
252  }
253  Time_.ResetStartTime();
254  // Loop through the vectors
255  for(int vecNum = 0; vecNum < NumVectors; vecNum++){
256  CallEuclid(X[vecNum], Y[vecNum]);
257  }
258  if(SetStats_ != 0){
259  Euclid_dhPrintTestData(eu, stdout);
260  }
261  NumApplyInverse_ = NumApplyInverse_ + 1;
262  ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
263  return 0;
264 } //ApplyInverse()
265 
266 //==============================================================================
267 int Ifpack_Euclid::CallEuclid(double *x, double *y) const{
268  Euclid_dhApply(eu, x, y);
269  return 0;
270 } //CallEuclid()
271 
272 //==============================================================================
273 std::ostream& operator << (std::ostream& os, const Ifpack_Euclid& A){
274  if (!A.Comm().MyPID()) {
275  os << endl;
276  os << "================================================================================" << endl;
277  os << "Ifpack_Euclid: " << A.Label () << endl << endl;
278  os << "Using " << A.Comm().NumProc() << " processors." << endl;
279  os << "Global number of rows = " << A.Matrix().NumGlobalRows() << endl;
280  os << "Global number of nonzeros = " << A.Matrix().NumGlobalNonzeros() << endl;
281  os << "Condition number estimate = " << A.Condest() << endl;
282  os << endl;
283  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
284  os << "----- ------- -------------- ------------ --------" << endl;
285  os << "Initialize() " << std::setw(5) << A.NumInitialize()
286  << " " << std::setw(15) << A.InitializeTime()
287  << " 0.0 0.0" << endl;
288  os << "Compute() " << std::setw(5) << A.NumCompute()
289  << " " << std::setw(15) << A.ComputeTime()
290  << " " << std::setw(15) << 1.0e-6 * A.ComputeFlops();
291  if (A.ComputeTime() != 0.0)
292  os << " " << std::setw(15) << 1.0e-6 * A.ComputeFlops() / A.ComputeTime() << endl;
293  else
294  os << " " << std::setw(15) << 0.0 << endl;
295  os << "ApplyInverse() " << std::setw(5) << A.NumApplyInverse()
296  << " " << std::setw(15) << A.ApplyInverseTime()
297  << " " << std::setw(15) << 1.0e-6 * A.ApplyInverseFlops();
298  if (A.ApplyInverseTime() != 0.0)
299  os << " " << std::setw(15) << 1.0e-6 * A.ApplyInverseFlops() / A.ApplyInverseTime() << endl;
300  else
301  os << " " << std::setw(15) << 0.0 << endl;
302  os << "================================================================================" << endl;
303  os << endl;
304  }
305  return os;
306 } // <<
307 
308 //==============================================================================
309 double Ifpack_Euclid::Condest(const Ifpack_CondestType CT,
310  const int MaxIters,
311  const double Tol,
312  Epetra_RowMatrix* Matrix_in){
313  if (!IsComputed()) // cannot compute right now
314  return(-1.0);
315  return(Condest_);
316 } //Condest() - not implemented
317 
318 #endif // HAVE_EUCLID && HAVE_MPI