IFPACK  Development
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 //==============================================================================
94 Ifpack_ILUT::~Ifpack_ILUT()
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  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 //==========================================================================
144 int Ifpack_ILUT::Initialize()
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>
180 int Ifpack_ILUT::TCompute()
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 
495 int Ifpack_ILUT::Compute() {
496  if (!IsInitialized())
497  IFPACK_CHK_ERR(Initialize());
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 
573  ++NumApplyInverse_;
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 //=============================================================================
591 double Ifpack_ILUT::Condest(const Ifpack_CondestType CT,
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 }
Ifpack_ILUT::ApplyInverse
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.
Definition: Ifpack_ILUT.cpp:537
Ifpack_HashTable
Definition: Ifpack_HashTable.h:258
Epetra_Time::ResetStartTime
void ResetStartTime(void)
Epetra_Comm::NumProc
virtual int NumProc() const=0
Epetra_MultiVector::Random
int Random()
Ifpack_ILUT::ApplyInverseTime
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition: Ifpack_ILUT.h:248
Ifpack_Utils.h
Epetra_RowMatrix::NumGlobalNonzeros
virtual int NumGlobalNonzeros() const=0
Ifpack_ILUT::Print
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
Definition: Ifpack_ILUT.cpp:607
Ifpack_ILUT::IsInitialized
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition: Ifpack_ILUT.h:113
Ifpack_AbsComp
Definition: Ifpack_IKLU.cpp:217
Ifpack_ILUT::NumInitialize
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack_ILUT.h:218
Ifpack_ILUT::ApplyInverseFlops
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition: Ifpack_ILUT.h:264
Ifpack_ILUT::IsComputed
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition: Ifpack_ILUT.h:130
Ifpack_ILUT::Condest
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
Definition: Ifpack_ILUT.h:154
Ifpack_ILUT::SetParameters
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
Definition: Ifpack_ILUT.cpp:107
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const=0
Epetra_CrsMatrix
Epetra_RowMatrix::NumMyRows
virtual int NumMyRows() const=0
Epetra_MultiVector::Pointers
double ** Pointers() const
Ifpack_ILUT::Label
const char * Label() const
Returns the label of this object.
Definition: Ifpack_ILUT.h:202
Ifpack_ILUT::Comm
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_ILUT.h:187
Epetra_RowMatrix
Epetra_SerialComm
Ifpack_ILUT::Initialize
int Initialize()
Initialize L and U with values from user matrix A.
Definition: Ifpack_ILUT.cpp:144
Ifpack_ILUT::DropTolerance
double DropTolerance() const
Gets the dropping tolerance.
Definition: Ifpack_ILUT.h:291
Epetra_Vector
Epetra_Comm::SumAll
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
Epetra_RowMatrix::ExtractMyRowCopy
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Ifpack_ILUT::RelativeThreshold
double RelativeThreshold() const
Get relative threshold value.
Definition: Ifpack_ILUT.h:285
Ifpack_ILUT::Compute
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
Definition: Ifpack_ILUT.cpp:495
Ifpack_ILUT::Ifpack_ILUT
Ifpack_ILUT(const Epetra_RowMatrix *A)
Ifpack_ILUT constuctor with variable number of indices per row.
Definition: Ifpack_ILUT.cpp:66
Ifpack_toString
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
Definition: Ifpack_Utils.cpp:242
Epetra_RowMatrix::MaxNumEntries
virtual int MaxNumEntries() const=0
Epetra_MultiVector
Ifpack_ILUT::InitializeTime
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack_ILUT.h:236
Ifpack_ILUT::ComputeFlops
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack_ILUT.h:259
Ifpack_ILUT::~Ifpack_ILUT
virtual ~Ifpack_ILUT()
Ifpack_ILUT Destructor.
Definition: Ifpack_ILUT.cpp:94
Epetra_Time::ElapsedTime
double ElapsedTime(void) const
Ifpack_ILUT::Matrix
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack_ILUT.h:190
Ifpack_ILUT::NumCompute
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack_ILUT.h:224
Ifpack_ILUT::NumApplyInverse
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition: Ifpack_ILUT.h:230
Ifpack_ILUT::AbsoluteThreshold
double AbsoluteThreshold() const
Get absolute threshold value.
Definition: Ifpack_ILUT.h:279
Ifpack_ILUT::RelaxValue
double RelaxValue() const
Set relative threshold value.
Definition: Ifpack_ILUT.h:274
Ifpack_ILUT::ComputeTime
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack_ILUT.h:242
Epetra_Map
Epetra_RowMatrix::Multiply
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
Epetra_MultiVector::NumVectors
int NumVectors() const
Generated by   doxygen 1.8.16