Ifpack Package Browser (Single Doxygen Collection)  Development
Ifpack_ILUT.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_Preconditioner.h"
45 #include "Ifpack_ILUT.h"
46 #include "Ifpack_Condest.h"
47 #include "Ifpack_Utils.h"
48 #include "Ifpack_HashTable.h"
49 #include "Epetra_SerialComm.h"
50 #include "Epetra_Comm.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_RowMatrix.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Epetra_Vector.h"
55 #include "Epetra_MultiVector.h"
56 #include "Epetra_Util.h"
57 #include "Teuchos_ParameterList.hpp"
58 #include "Teuchos_RefCountPtr.hpp"
59 #include <functional>
60 #include <algorithm>
61 
62 using namespace Teuchos;
63 
64 //==============================================================================
65 // FIXME: allocate Comm_ and Time_ the first Initialize() call
66 Ifpack_ILUT::Ifpack_ILUT(const Epetra_RowMatrix* A) :
67  A_(*A),
68  Comm_(A->Comm()),
69  Condest_(-1.0),
70  Relax_(0.),
71  Athresh_(0.0),
72  Rthresh_(1.0),
73  LevelOfFill_(1.0),
74  DropTolerance_(1e-12),
75  IsInitialized_(false),
76  IsComputed_(false),
77  UseTranspose_(false),
78  NumMyRows_(-1),
79  NumInitialize_(0),
80  NumCompute_(0),
81  NumApplyInverse_(0),
82  InitializeTime_(0.0),
83  ComputeTime_(0.0),
84  ApplyInverseTime_(0.0),
85  ComputeFlops_(0.0),
86  ApplyInverseFlops_(0.0),
87  Time_(Comm()),
88  GlobalNonzeros_(0)
89 {
90  // do nothing here..
91 }
92 
93 //==============================================================================
95 {
96  Destroy();
97 }
98 
99 //==============================================================================
101 {
102  IsInitialized_ = false;
103  IsComputed_ = false;
104 }
105 
106 //==========================================================================
108 {
109  using std::cerr;
110  using std::endl;
111 
112  try
113  {
114  LevelOfFill_ = List.get<double>("fact: ilut level-of-fill", LevelOfFill());
115  if (LevelOfFill_ <= 0.0)
116  IFPACK_CHK_ERR(-2); // must be greater than 0.0
117 
118  Athresh_ = List.get<double>("fact: absolute threshold", Athresh_);
119  Rthresh_ = List.get<double>("fact: relative threshold", Rthresh_);
120  Relax_ = List.get<double>("fact: relax value", Relax_);
121  DropTolerance_ = List.get<double>("fact: drop tolerance", DropTolerance_);
122 
123  Label_ = "IFPACK ILUT (fill=" + Ifpack_toString(LevelOfFill())
124  + ", relax=" + Ifpack_toString(RelaxValue())
125  + ", athr=" + Ifpack_toString(AbsoluteThreshold())
126  + ", rthr=" + Ifpack_toString(RelativeThreshold())
127  + ", droptol=" + Ifpack_toString(DropTolerance())
128  + ")";
129  return(0);
130  }
131  catch (...)
132  {
133  cerr << "Caught an exception while parsing the parameter list" << endl;
134  cerr << "This typically means that a parameter was set with the" << endl;
135  cerr << "wrong type (for example, int instead of double). " << endl;
136  cerr << "please check the documentation for the type required by each parameter." << endl;
137  IFPACK_CHK_ERR(-1);
138  }
139 
140  //return(0); // unreachable
141 }
142 
143 //==========================================================================
145 {
146  // delete previously allocated factorization
147  Destroy();
148 
149  Time_.ResetStartTime();
150 
151  // check only in serial
152  if (Comm().NumProc() == 1 && Matrix().NumMyRows() != Matrix().NumMyCols())
153  IFPACK_CHK_ERR(-2);
154 
155  NumMyRows_ = Matrix().NumMyRows();
156 
157  // nothing else to do here
158  IsInitialized_ = true;
159  ++NumInitialize_;
160  InitializeTime_ += Time_.ElapsedTime();
161 
162  return(0);
163 }
164 
165 //==========================================================================
166 class Ifpack_AbsComp
167 {
168  public:
169  inline bool operator()(const double& x, const double& y)
170  {
171  return(IFPACK_ABS(x) > IFPACK_ABS(y));
172  }
173 };
174 
175 //==========================================================================
176 // MS // completely rewritten the algorithm on 25-Jan-05, using more STL
177 // MS // and in a nicer and cleaner way. Also, it is more efficient.
178 // MS // WARNING: Still not fully tested!
179 template<typename int_type>
181 {
182  int Length = A_.MaxNumEntries();
183  std::vector<double> RowValuesL(Length);
184  std::vector<int> RowIndicesU(Length);
185  std::vector<int_type> RowIndicesU_LL;
186  RowIndicesU_LL.resize(Length);
187  std::vector<double> RowValuesU(Length);
188 
189  int RowNnzU;
190 
191  L_ = rcp(new Epetra_CrsMatrix(Copy, *SerialMap_, 0));
192  U_ = rcp(new Epetra_CrsMatrix(Copy, *SerialMap_, 0));
193 
194  if ((L_.get() == 0) || (U_.get() == 0))
195  IFPACK_CHK_ERR(-5); // memory allocation error
196 
197  // insert first row in U_ and L_
198  IFPACK_CHK_ERR(A_.ExtractMyRowCopy(0,Length,RowNnzU,
199  &RowValuesU[0],&RowIndicesU[0]));
200 
201  bool distributed = (Comm().NumProc() > 1)?true:false;
202 
203  if (distributed)
204  {
205  int count = 0;
206  for (int i = 0 ;i < RowNnzU ; ++i)
207  {
208  if (RowIndicesU[i] < NumMyRows_){
209  RowIndicesU[count] = RowIndicesU[i];
210  RowValuesU[count] = RowValuesU[i];
211  ++count;
212  }
213  else
214  continue;
215  }
216  RowNnzU = count;
217  }
218 
219  // modify diagonal
220  for (int i = 0 ;i < RowNnzU ; ++i) {
221  if (RowIndicesU[i] == 0) {
222  double& v = RowValuesU[i];
223  v = AbsoluteThreshold() * EPETRA_SGN(v) + RelativeThreshold() * v;
224  break;
225  }
226  }
227 
228  std::copy(&(RowIndicesU[0]), &(RowIndicesU[0]) + RowNnzU, RowIndicesU_LL.begin());
229  IFPACK_CHK_ERR(U_->InsertGlobalValues(0,RowNnzU,&(RowValuesU[0]),
230  &(RowIndicesU_LL[0])));
231  // FIXME: DOES IT WORK IN PARALLEL ??
232  RowValuesU[0] = 1.0;
233  RowIndicesU[0] = 0;
234  int_type RowIndicesU_0_templ = RowIndicesU[0];
235  IFPACK_CHK_ERR(L_->InsertGlobalValues(0,1,&(RowValuesU[0]),
236  &RowIndicesU_0_templ));
237 
238  int max_keys = (int) 10 * A_.MaxNumEntries() * LevelOfFill() ;
239  Ifpack_HashTable SingleRowU(max_keys , 1);
240  Ifpack_HashTable SingleRowL(max_keys , 1);
241 
242  int hash_size = SingleRowU.getRecommendedHashSize(max_keys) ;
243  std::vector<int> keys; keys.reserve(hash_size * 10);
244  std::vector<double> values; values.reserve(hash_size * 10);
245  std::vector<double> AbsRow; AbsRow.reserve(hash_size * 10);
246 
247  // =================== //
248  // start factorization //
249  // =================== //
250 
251 #ifdef IFPACK_FLOPCOUNTERS
252  double this_proc_flops = 0.0;
253 #endif
254 
255  for (int row_i = 1 ; row_i < NumMyRows_ ; ++row_i)
256  {
257  // get row `row_i' of the matrix, store in U pointers
258  IFPACK_CHK_ERR(A_.ExtractMyRowCopy(row_i,Length,RowNnzU,
259  &RowValuesU[0],&RowIndicesU[0]));
260 
261  if (distributed)
262  {
263  int count = 0;
264  for (int i = 0 ;i < RowNnzU ; ++i)
265  {
266  if (RowIndicesU[i] < NumMyRows_){
267  RowIndicesU[count] = RowIndicesU[i];
268  RowValuesU[count] = RowValuesU[i];
269  ++count;
270  }
271  else
272  continue;
273  }
274  RowNnzU = count;
275  }
276 
277  int NnzLower = 0;
278  int NnzUpper = 0;
279 
280  for (int i = 0 ;i < RowNnzU ; ++i) {
281  if (RowIndicesU[i] < row_i)
282  NnzLower++;
283  else if (RowIndicesU[i] == row_i) {
284  // add threshold
285  NnzUpper++;
286  double& v = RowValuesU[i];
287  v = AbsoluteThreshold() * EPETRA_SGN(v) + RelativeThreshold() * v;
288  }
289  else
290  NnzUpper++;
291  }
292 
293  int FillL = (int)(LevelOfFill() * NnzLower);
294  int FillU = (int)(LevelOfFill() * NnzUpper);
295  if (FillL == 0) FillL = 1;
296  if (FillU == 0) FillU = 1;
297 
298  // convert line `row_i' into STL map for fast access
299  SingleRowU.reset();
300 
301  for (int i = 0 ; i < RowNnzU ; ++i) {
302  SingleRowU.set(RowIndicesU[i], RowValuesU[i]);
303  }
304 
305  // for the multipliers
306  SingleRowL.reset();
307 
308  int start_col = NumMyRows_;
309  for (int i = 0 ; i < RowNnzU ; ++i)
310  start_col = EPETRA_MIN(start_col, RowIndicesU[i]);
311 
312 #ifdef IFPACK_FLOPCOUNTERS
313  int flops = 0;
314 #endif
315 
316  for (int col_k = start_col ; col_k < row_i ; ++col_k) {
317 
318  int_type* ColIndicesK;
319  double* ColValuesK;
320  int ColNnzK;
321 
322  double xxx = SingleRowU.get(col_k);
323  // This factorization is too "relaxed". Leaving it as it is, as Tifpack
324  // does it correctly.
325  if (IFPACK_ABS(xxx) > DropTolerance()) {
326  IFPACK_CHK_ERR(U_->ExtractGlobalRowView(col_k, ColNnzK, ColValuesK,
327  ColIndicesK));
328 
329  // FIXME: can keep trace of diagonals
330  double DiagonalValueK = 0.0;
331  for (int i = 0 ; i < ColNnzK ; ++i) {
332  if (ColIndicesK[i] == col_k) {
333  DiagonalValueK = ColValuesK[i];
334  break;
335  }
336  }
337 
338  // FIXME: this should never happen!
339  if (DiagonalValueK == 0.0)
340  DiagonalValueK = AbsoluteThreshold();
341 
342  double yyy = xxx / DiagonalValueK ;
343  SingleRowL.set(col_k, yyy);
344 #ifdef IFPACK_FLOPCOUNTERS
345  ++flops;
346 #endif
347 
348  for (int j = 0 ; yyy != 0.0 && j < ColNnzK ; ++j)
349  {
350  int_type col_j = ColIndicesK[j];
351 
352  if (col_j < col_k) continue;
353 
354  SingleRowU.set(col_j, -yyy * ColValuesK[j], true);
355 #ifdef IFPACK_FLOPCOUNTERS
356  flops += 2;
357 #endif
358  }
359  }
360  }
361 
362 #ifdef IFPACK_FLOPCOUNTERS
363  this_proc_flops += 1.0 * flops;
364 #endif
365 
366  double cutoff = DropTolerance();
367  double DiscardedElements = 0.0;
368  int count;
369 
370  // drop elements to satisfy LevelOfFill(), start with L
371  count = 0;
372  int sizeL = SingleRowL.getNumEntries();
373  keys.resize(sizeL);
374  values.resize(sizeL);
375 
376  AbsRow.resize(sizeL);
377 
378  SingleRowL.arrayify(
379  keys.size() ? &keys[0] : 0,
380  values.size() ? &values[0] : 0
381  );
382  for (int i = 0; i < sizeL; ++i)
383  if (IFPACK_ABS(values[i]) > DropTolerance()) {
384  AbsRow[count++] = IFPACK_ABS(values[i]);
385  }
386 
387  if (count > FillL) {
388  nth_element(AbsRow.begin(), AbsRow.begin() + FillL, AbsRow.begin() + count,
389  std::greater<double>());
390  cutoff = AbsRow[FillL];
391  }
392 
393  for (int i = 0; i < sizeL; ++i) {
394  if (IFPACK_ABS(values[i]) >= cutoff) {
395  int_type key_templ = keys[i];
396  IFPACK_CHK_ERR(L_->InsertGlobalValues(row_i,1, &values[i], &key_templ));
397  }
398  else
399  DiscardedElements += values[i];
400  }
401 
402  // FIXME: DOES IT WORK IN PARALLEL ???
403  // add 1 to the diagonal
404  double dtmp = 1.0;
405  int_type row_i_templ = row_i;
406  IFPACK_CHK_ERR(L_->InsertGlobalValues(row_i,1, &dtmp, &row_i_templ));
407 
408  // same business with U_
409  count = 0;
410  int sizeU = SingleRowU.getNumEntries();
411  AbsRow.resize(sizeU + 1);
412  keys.resize(sizeU + 1);
413  values.resize(sizeU + 1);
414 
415  SingleRowU.arrayify(&keys[0], &values[0]);
416 
417  for (int i = 0; i < sizeU; ++i)
418  if (keys[i] >= row_i && IFPACK_ABS(values[i]) > DropTolerance())
419  {
420  AbsRow[count++] = IFPACK_ABS(values[i]);
421  }
422 
423  if (count > FillU) {
424  nth_element(AbsRow.begin(), AbsRow.begin() + FillU, AbsRow.begin() + count,
425  std::greater<double>());
426  cutoff = AbsRow[FillU];
427  }
428 
429  // sets the factors in U_
430  for (int i = 0; i < sizeU; ++i)
431  {
432  int col = keys[i];
433  double val = values[i];
434 
435  if (col >= row_i) {
436  if (IFPACK_ABS(val) >= cutoff || row_i == col) {
437  int_type col_templ = col;
438  IFPACK_CHK_ERR(U_->InsertGlobalValues(row_i,1, &val, &col_templ));
439  }
440  else
441  DiscardedElements += val;
442  }
443  }
444 
445  // FIXME: not so sure of that!
446  if (RelaxValue() != 0.0) {
447  DiscardedElements *= RelaxValue();
448  IFPACK_CHK_ERR(U_->InsertGlobalValues(row_i,1, &DiscardedElements,
449  &row_i_templ));
450  }
451  }
452 
453 #ifdef IFPACK_FLOPCOUNTERS
454  double tf;
455  Comm().SumAll(&this_proc_flops, &tf, 1);
456  ComputeFlops_ += tf;
457 #endif
458 
459  IFPACK_CHK_ERR(L_->FillComplete());
460  IFPACK_CHK_ERR(U_->FillComplete());
461 
462 #if 0
463  // to check the complete factorization
464  Epetra_Vector LHS(A_.RowMatrixRowMap());
465  Epetra_Vector RHS1(A_.RowMatrixRowMap());
466  Epetra_Vector RHS2(A_.RowMatrixRowMap());
467  Epetra_Vector RHS3(A_.RowMatrixRowMap());
468  LHS.Random();
469 
470  cout << "A = " << A_.NumGlobalNonzeros() << endl;
471  cout << "L = " << L_->NumGlobalNonzeros() << endl;
472  cout << "U = " << U_->NumGlobalNonzeros() << endl;
473 
474  A_.Multiply(false,LHS,RHS1);
475  U_->Multiply(false,LHS,RHS2);
476  L_->Multiply(false,RHS2,RHS3);
477 
478  RHS1.Update(-1.0, RHS3, 1.0);
479  double Norm;
480  RHS1.Norm2(&Norm);
481 #endif
482 
483  long long MyNonzeros = L_->NumGlobalNonzeros64() + U_->NumGlobalNonzeros64();
484  Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
485 
486  IsComputed_ = true;
487 
488  ++NumCompute_;
489  ComputeTime_ += Time_.ElapsedTime();
490 
491  return(0);
492 
493 }
494 
496  if (!IsInitialized())
498 
499  Time_.ResetStartTime();
500  IsComputed_ = false;
501 
502  NumMyRows_ = A_.NumMyRows();
503  bool distributed = (Comm().NumProc() > 1)?true:false;
504 
505 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
506  if (distributed)
507  {
508  SerialComm_ = rcp(new Epetra_SerialComm);
509  SerialMap_ = rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
510  assert (SerialComm_.get() != 0);
511  assert (SerialMap_.get() != 0);
512  }
513  else
514  SerialMap_ = rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);
515 #endif
516 
517  int ret_val = 1;
518 
519 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
520  if(SerialMap_->GlobalIndicesInt()) {
521  ret_val = TCompute<int>();
522  }
523  else
524 #endif
525 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
526  if(SerialMap_->GlobalIndicesLongLong()) {
527  ret_val = TCompute<long long>();
528  }
529  else
530 #endif
531  throw "Ifpack_ILUT::Compute: GlobalIndices type unknown";
532 
533  return ret_val;
534 }
535 
536 //=============================================================================
537 int Ifpack_ILUT::ApplyInverse(const Epetra_MultiVector& X,
538  Epetra_MultiVector& Y) const
539 {
540  if (!IsComputed())
541  IFPACK_CHK_ERR(-2); // compute preconditioner first
542 
543  if (X.NumVectors() != Y.NumVectors())
544  IFPACK_CHK_ERR(-3); // Return error: X and Y not the same size
545 
546  Time_.ResetStartTime();
547 
548  // NOTE: L_ and U_ are based on SerialMap_, while Xcopy is based
549  // on A.Map()... which are in general different. However, Solve()
550  // does not seem to care... which is fine with me.
551  //
552  // AztecOO gives X and Y pointing to the same memory location,
553  // need to create an auxiliary vector, Xcopy
554  Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
555  if (X.Pointers()[0] == Y.Pointers()[0])
556  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
557  else
558  Xcopy = Teuchos::rcp( &X, false );
559 
560  if (!UseTranspose_)
561  {
562  // solves LU Y = X
563  IFPACK_CHK_ERR(L_->Solve(false,false,false,*Xcopy,Y));
564  IFPACK_CHK_ERR(U_->Solve(true,false,false,Y,Y));
565  }
566  else
567  {
568  // solves U(trans) L(trans) Y = X
569  IFPACK_CHK_ERR(U_->Solve(true,true,false,*Xcopy,Y));
570  IFPACK_CHK_ERR(L_->Solve(false,true,false,Y,Y));
571  }
572 
574 #ifdef IFPACK_FLOPCOUNTERS
575  ApplyInverseFlops_ += X.NumVectors() * 2 * GlobalNonzeros_;
576 #endif
577  ApplyInverseTime_ += Time_.ElapsedTime();
578 
579  return(0);
580 
581 }
582 //=============================================================================
583 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
584 int Ifpack_ILUT::Apply(const Epetra_MultiVector& X,
585  Epetra_MultiVector& Y) const
586 {
587  return(-98);
588 }
589 
590 //=============================================================================
592  const int MaxIters, const double Tol,
593  Epetra_RowMatrix* Matrix_in)
594 {
595  if (!IsComputed()) // cannot compute right now
596  return(-1.0);
597 
598  // NOTE: this is computing the *local* condest
599  if (Condest_ == -1.0)
600  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
601 
602  return(Condest_);
603 }
604 
605 //=============================================================================
606 std::ostream&
607 Ifpack_ILUT::Print(std::ostream& os) const
608 {
609  using std::endl;
610 
611  if (!Comm().MyPID()) {
612  os << endl;
613  os << "================================================================================" << endl;
614  os << "Ifpack_ILUT: " << Label() << endl << endl;
615  os << "Level-of-fill = " << LevelOfFill() << endl;
616  os << "Absolute threshold = " << AbsoluteThreshold() << endl;
617  os << "Relative threshold = " << RelativeThreshold() << endl;
618  os << "Relax value = " << RelaxValue() << endl;
619  os << "Condition number estimate = " << Condest() << endl;
620  os << "Global number of rows = " << A_.NumGlobalRows64() << endl;
621  if (IsComputed_) {
622  os << "Number of nonzeros in A = " << A_.NumGlobalNonzeros64() << endl;
623  os << "Number of nonzeros in L + U = " << NumGlobalNonzeros64()
624  << " ( = " << 100.0 * NumGlobalNonzeros64() / A_.NumGlobalNonzeros64()
625  << " % of A)" << endl;
626  os << "nonzeros / rows = "
627  << 1.0 * NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
628  }
629  os << endl;
630  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
631  os << "----- ------- -------------- ------------ --------" << endl;
632  os << "Initialize() " << std::setw(5) << NumInitialize()
633  << " " << std::setw(15) << InitializeTime()
634  << " 0.0 0.0" << endl;
635  os << "Compute() " << std::setw(5) << NumCompute()
636  << " " << std::setw(15) << ComputeTime()
637  << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
638  if (ComputeTime() != 0.0)
639  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
640  else
641  os << " " << std::setw(15) << 0.0 << endl;
642  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
643  << " " << std::setw(15) << ApplyInverseTime()
644  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
645  if (ApplyInverseTime() != 0.0)
646  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
647  else
648  os << " " << std::setw(15) << 0.0 << endl;
649  os << "================================================================================" << endl;
650  os << endl;
651  }
652 
653  return(os);
654 }
double AbsoluteThreshold() const
Get absolute threshold value.
Definition: Ifpack_ILUT.h:279
double RelaxValue() const
Set relative threshold value.
Definition: Ifpack_ILUT.h:274
double DropTolerance() const
Gets the dropping tolerance.
Definition: Ifpack_ILUT.h:291
double get(const key_type key)
Returns an element from the hash table, or 0.0 if not found.
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
U factor.
Definition: Ifpack_ILUT.h:347
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack_ILUT.h:218
double RelativeThreshold() const
Get relative threshold value.
Definition: Ifpack_ILUT.h:285
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition: Ifpack_ILUT.h:113
T & get(ParameterList &l, const std::string &name)
double DropTolerance_
Discards all elements below this tolerance.
Definition: Ifpack_ILUT.h:359
Teuchos::RefCountPtr< Epetra_Map > SerialMap_
Definition: Ifpack_ILUT.h:391
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
L factor.
Definition: Ifpack_ILUT.h:345
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
bool IsInitialized_
true if this object has been initialized
Definition: Ifpack_ILUT.h:363
Ifpack_ILUT(const Epetra_RowMatrix *A)
Ifpack_ILUT constuctor with variable number of indices per row.
Definition: Ifpack_ILUT.cpp:66
long long GlobalNonzeros_
Global number of nonzeros in L and U factors.
Definition: Ifpack_ILUT.h:389
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition: Ifpack_ILUT.h:248
long long NumGlobalNonzeros64() const
Definition: Ifpack_ILUT.h:303
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
int Initialize()
Initialize L and U with values from user matrix A.
virtual ~Ifpack_ILUT()
Ifpack_ILUT Destructor.
Definition: Ifpack_ILUT.cpp:94
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition: Ifpack_ILUT.h:130
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
bool UseTranspose_
true if transpose has to be used.
Definition: Ifpack_ILUT.h:367
double Condest_
Condition number estimate.
Definition: Ifpack_ILUT.h:349
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int NumInitialize_
Contains the number of successful calls to Initialize().
Definition: Ifpack_ILUT.h:371
std::string Label_
Label for this object.
Definition: Ifpack_ILUT.h:361
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
int getNumEntries() const
Returns the number of stored entries.
void set(const key_type key, const double value, const bool addToValue=false)
Sets an element in the hash table.
Teuchos::RefCountPtr< Epetra_SerialComm > SerialComm_
Definition: Ifpack_ILUT.h:390
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack_ILUT.h:190
double ComputeTime_
Contains the time for all successful calls to Compute().
Definition: Ifpack_ILUT.h:379
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack_ILUT.h:236
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition: Ifpack_ILUT.h:264
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
#define false
const Epetra_RowMatrix & A_
reference to the matrix to be preconditioned.
Definition: Ifpack_ILUT.h:341
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack_ILUT.h:224
bool IsComputed_
true if this object has been computed
Definition: Ifpack_ILUT.h:365
void arrayify(key_type *key_array, double *val_array)
Converts the contents in array format for both keys and values.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Definition: Ifpack_ILUT.h:381
double LevelOfFill_
Level-of-fill.
Definition: Ifpack_ILUT.h:357
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition: Ifpack_ILUT.h:230
#define IFPACK_ABS(x)
const char * Label() const
Returns the label of this object.
Definition: Ifpack_ILUT.h:202
void Destroy()
Releases all allocated memory.
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
Definition: Ifpack_ILUT.h:154
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_ILUT forward/back solve on a Epetra_MultiVector X in Y...
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Definition: Ifpack_ILUT.h:375
double ComputeFlops_
Contains the number of flops for Compute().
Definition: Ifpack_ILUT.h:383
Copy
double LevelOfFill() const
Definition: Ifpack_ILUT.h:269
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
Definition: Ifpack_ILUT.h:385
double Rthresh_
Relative threshold.
Definition: Ifpack_ILUT.h:355
bool operator()(const double &x, const double &y)
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_ILUT.h:187
double InitializeTime_
Contains the time for all successful calls to Initialize().
Definition: Ifpack_ILUT.h:377
int NumCompute_
Contains the number of successful call to Compute().
Definition: Ifpack_ILUT.h:373
Epetra_Time Time_
Used for timing purposed.
Definition: Ifpack_ILUT.h:387
int NumMyRows_
Number of local rows.
Definition: Ifpack_ILUT.h:369
#define IFPACK_CHK_ERR(ifpack_err)
int getRecommendedHashSize(int n)
double Athresh_
Absolute threshold.
Definition: Ifpack_ILUT.h:353
void reset()
Resets the entries of the already allocated memory. This method can be used to clean an array...
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack_ILUT.h:259
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack_ILUT.h:242
double Relax_
relaxation value
Definition: Ifpack_ILUT.h:351