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