IFPACK Development
Loading...
Searching...
No Matches
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
62using namespace Teuchos;
63
64//==============================================================================
65// FIXME: allocate Comm_ and Time_ the first Initialize() call
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//==============================================================================
100void Ifpack_ILUT::Destroy()
101{
102 IsInitialized_ = false;
103 IsComputed_ = false;
104}
105
106//==========================================================================
107int 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//==========================================================================
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//==========================================================================
166class 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!
179template<typename int_type>
180int 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
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())
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//=============================================================================
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
584int Ifpack_ILUT::Apply(const Epetra_MultiVector& /* X */,
585 Epetra_MultiVector& /* Y */) const
586{
587 return(-98);
588}
589
590//=============================================================================
591double 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//=============================================================================
606std::ostream&
607Ifpack_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}
virtual int NumProc() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int NumVectors() const
double ** Pointers() const
virtual int NumMyRows() const=0
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual int NumGlobalNonzeros() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
void ResetStartTime(void)
double ElapsedTime(void) const
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
double RelativeThreshold() const
Get relative threshold value.
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
virtual double InitializeTime() const
Returns the time spent in Initialize().
double DropTolerance() const
Gets the dropping tolerance.
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 RelaxValue() const
Set relative threshold value.
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Ifpack_ILUT(const Epetra_RowMatrix *A)
Ifpack_ILUT constuctor with variable number of indices per row.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
const char * Label() const
Returns the label of this object.
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
int Initialize()
Initialize L and U with values from user matrix A.
virtual double ComputeTime() const
Returns the time spent in Compute().
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual ~Ifpack_ILUT()
Ifpack_ILUT Destructor.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
double AbsoluteThreshold() const
Get absolute threshold value.