43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_Preconditioner.h"
45 #include "Ifpack_ILUT.h"
46 #include "Ifpack_Condest.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"
62 using namespace Teuchos;
74 DropTolerance_(1e-12),
75 IsInitialized_(false),
84 ApplyInverseTime_(0.0),
86 ApplyInverseFlops_(0.0),
100 void Ifpack_ILUT::Destroy()
102 IsInitialized_ =
false;
114 LevelOfFill_ = List.get<
double>(
"fact: ilut level-of-fill", LevelOfFill());
115 if (LevelOfFill_ <= 0.0)
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_);
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;
158 IsInitialized_ =
true;
169 inline bool operator()(
const double& x,
const double& y)
171 return(IFPACK_ABS(x) > IFPACK_ABS(y));
179 template<
typename int_type>
180 int Ifpack_ILUT::TCompute()
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);
194 if ((L_.get() == 0) || (U_.get() == 0))
199 &RowValuesU[0],&RowIndicesU[0]));
201 bool distributed = (
Comm().
NumProc() > 1)?
true:
false;
206 for (
int i = 0 ;i < RowNnzU ; ++i)
208 if (RowIndicesU[i] < NumMyRows_){
209 RowIndicesU[count] = RowIndicesU[i];
210 RowValuesU[count] = RowValuesU[i];
220 for (
int i = 0 ;i < RowNnzU ; ++i) {
221 if (RowIndicesU[i] == 0) {
222 double& v = RowValuesU[i];
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])));
234 int_type RowIndicesU_0_templ = RowIndicesU[0];
235 IFPACK_CHK_ERR(L_->InsertGlobalValues(0,1,&(RowValuesU[0]),
236 &RowIndicesU_0_templ));
238 int max_keys = (int) 10 * A_.
MaxNumEntries() * LevelOfFill() ;
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);
251 #ifdef IFPACK_FLOPCOUNTERS
252 double this_proc_flops = 0.0;
255 for (
int row_i = 1 ; row_i < NumMyRows_ ; ++row_i)
259 &RowValuesU[0],&RowIndicesU[0]));
264 for (
int i = 0 ;i < RowNnzU ; ++i)
266 if (RowIndicesU[i] < NumMyRows_){
267 RowIndicesU[count] = RowIndicesU[i];
268 RowValuesU[count] = RowValuesU[i];
280 for (
int i = 0 ;i < RowNnzU ; ++i) {
281 if (RowIndicesU[i] < row_i)
283 else if (RowIndicesU[i] == row_i) {
286 double& v = RowValuesU[i];
293 int FillL = (int)(LevelOfFill() * NnzLower);
294 int FillU = (int)(LevelOfFill() * NnzUpper);
295 if (FillL == 0) FillL = 1;
296 if (FillU == 0) FillU = 1;
301 for (
int i = 0 ; i < RowNnzU ; ++i) {
302 SingleRowU.set(RowIndicesU[i], RowValuesU[i]);
308 int start_col = NumMyRows_;
309 for (
int i = 0 ; i < RowNnzU ; ++i)
310 start_col = EPETRA_MIN(start_col, RowIndicesU[i]);
312 #ifdef IFPACK_FLOPCOUNTERS
316 for (
int col_k = start_col ; col_k < row_i ; ++col_k) {
318 int_type* ColIndicesK;
322 double xxx = SingleRowU.get(col_k);
326 IFPACK_CHK_ERR(U_->ExtractGlobalRowView(col_k, ColNnzK, ColValuesK,
330 double DiagonalValueK = 0.0;
331 for (
int i = 0 ; i < ColNnzK ; ++i) {
332 if (ColIndicesK[i] == col_k) {
333 DiagonalValueK = ColValuesK[i];
339 if (DiagonalValueK == 0.0)
342 double yyy = xxx / DiagonalValueK ;
343 SingleRowL.set(col_k, yyy);
344 #ifdef IFPACK_FLOPCOUNTERS
348 for (
int j = 0 ; yyy != 0.0 && j < ColNnzK ; ++j)
350 int_type col_j = ColIndicesK[j];
352 if (col_j < col_k)
continue;
354 SingleRowU.set(col_j, -yyy * ColValuesK[j],
true);
355 #ifdef IFPACK_FLOPCOUNTERS
362 #ifdef IFPACK_FLOPCOUNTERS
363 this_proc_flops += 1.0 * flops;
367 double DiscardedElements = 0.0;
372 int sizeL = SingleRowL.getNumEntries();
374 values.resize(sizeL);
376 AbsRow.resize(sizeL);
379 keys.size() ? &keys[0] : 0,
380 values.size() ? &values[0] : 0
382 for (
int i = 0; i < sizeL; ++i)
384 AbsRow[count++] = IFPACK_ABS(values[i]);
388 nth_element(AbsRow.begin(), AbsRow.begin() + FillL, AbsRow.begin() + count,
389 std::greater<double>());
390 cutoff = AbsRow[FillL];
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));
399 DiscardedElements += values[i];
405 int_type row_i_templ = row_i;
406 IFPACK_CHK_ERR(L_->InsertGlobalValues(row_i,1, &dtmp, &row_i_templ));
410 int sizeU = SingleRowU.getNumEntries();
411 AbsRow.resize(sizeU + 1);
412 keys.resize(sizeU + 1);
413 values.resize(sizeU + 1);
415 SingleRowU.arrayify(&keys[0], &values[0]);
417 for (
int i = 0; i < sizeU; ++i)
418 if (keys[i] >= row_i && IFPACK_ABS(values[i]) >
DropTolerance())
420 AbsRow[count++] = IFPACK_ABS(values[i]);
424 nth_element(AbsRow.begin(), AbsRow.begin() + FillU, AbsRow.begin() + count,
425 std::greater<double>());
426 cutoff = AbsRow[FillU];
430 for (
int i = 0; i < sizeU; ++i)
433 double val = values[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));
441 DiscardedElements += val;
448 IFPACK_CHK_ERR(U_->InsertGlobalValues(row_i,1, &DiscardedElements,
453 #ifdef IFPACK_FLOPCOUNTERS
459 IFPACK_CHK_ERR(L_->FillComplete());
460 IFPACK_CHK_ERR(U_->FillComplete());
471 cout <<
"L = " << L_->NumGlobalNonzeros() << endl;
472 cout <<
"U = " << U_->NumGlobalNonzeros() << endl;
475 U_->Multiply(
false,LHS,RHS2);
476 L_->Multiply(
false,RHS2,RHS3);
478 RHS1.Update(-1.0, RHS3, 1.0);
483 long long MyNonzeros = L_->NumGlobalNonzeros64() + U_->NumGlobalNonzeros64();
484 Comm().
SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
503 bool distributed = (
Comm().
NumProc() > 1)?
true:
false;
505 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
509 SerialMap_ = rcp(
new Epetra_Map(NumMyRows_, 0, *SerialComm_));
510 assert (SerialComm_.get() != 0);
511 assert (SerialMap_.get() != 0);
514 SerialMap_ = rcp(const_cast<Epetra_Map*>(&A_.
RowMatrixRowMap()),
false);
519 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
520 if(SerialMap_->GlobalIndicesInt()) {
521 ret_val = TCompute<int>();
525 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
526 if(SerialMap_->GlobalIndicesLongLong()) {
527 ret_val = TCompute<long long>();
531 throw "Ifpack_ILUT::Compute: GlobalIndices type unknown";
543 if (X.NumVectors() != Y.NumVectors())
554 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
555 if (X.Pointers()[0] == Y.Pointers()[0])
558 Xcopy = Teuchos::rcp( &X,
false );
563 IFPACK_CHK_ERR(L_->Solve(
false,
false,
false,*Xcopy,Y));
564 IFPACK_CHK_ERR(U_->Solve(
true,
false,
false,Y,Y));
569 IFPACK_CHK_ERR(U_->Solve(
true,
true,
false,*Xcopy,Y));
570 IFPACK_CHK_ERR(L_->Solve(
false,
true,
false,Y,Y));
574 #ifdef IFPACK_FLOPCOUNTERS
575 ApplyInverseFlops_ += X.NumVectors() * 2 * GlobalNonzeros_;
592 const int MaxIters,
const double Tol,
599 if (Condest_ == -1.0)
600 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
611 if (!
Comm().MyPID()) {
613 os <<
"================================================================================" << endl;
614 os <<
"Ifpack_ILUT: " <<
Label() << endl << endl;
615 os <<
"Level-of-fill = " << LevelOfFill() << endl;
618 os <<
"Relax value = " <<
RelaxValue() << endl;
619 os <<
"Condition number estimate = " <<
Condest() << endl;
620 os <<
"Global number of rows = " << A_.NumGlobalRows64() << endl;
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;
630 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
631 os <<
"----- ------- -------------- ------------ --------" << endl;
634 <<
" 0.0 0.0" << endl;
635 os <<
"Compute() " << std::setw(5) <<
NumCompute()
641 os <<
" " << std::setw(15) << 0.0 << endl;
648 os <<
" " << std::setw(15) << 0.0 << endl;
649 os <<
"================================================================================" << endl;
virtual std::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.
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
virtual const Epetra_Map & RowMatrixRowMap() const =0
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
double ElapsedTime(void) const
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.
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
int Initialize()
Initialize L and U with values from user matrix A.
virtual int NumGlobalNonzeros() const =0
virtual ~Ifpack_ILUT()
Ifpack_ILUT Destructor.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
virtual int MaxNumEntries() const =0
double RelativeThreshold() const
Get relative threshold value.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual int NumCompute() const
Returns the number of calls to Compute().
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
virtual int NumMyRows() const =0
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
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...
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
double DropTolerance() const
Gets the dropping tolerance.
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual int NumProc() const =0
double AbsoluteThreshold() const
Get absolute threshold value.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
const char * Label() const
Returns the label of this object.
virtual double InitializeTime() const
Returns the time spent in Initialize().
void ResetStartTime(void)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
double RelaxValue() const
Set relative threshold value.