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"
74 DropTolerance_(1e-12),
75 IsInitialized_(
false),
84 ApplyInverseTime_(0.0),
86 ApplyInverseFlops_(0.0),
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;
179template<
typename int_type>
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)
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());
230 &(RowIndicesU_LL[0])));
234 int_type RowIndicesU_0_templ = RowIndicesU[0];
236 &RowIndicesU_0_templ));
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)
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];
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]);
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);
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;
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)
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) {
395 int_type key_templ = keys[i];
399 DiscardedElements += values[i];
405 int_type row_i_templ = row_i;
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)
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;
441 DiscardedElements += val;
453#ifdef IFPACK_FLOPCOUNTERS
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();
503 bool distributed = (
Comm().
NumProc() > 1)?
true:
false;
505#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
519#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
521 ret_val = TCompute<int>();
525#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
527 ret_val = TCompute<long long>();
531 throw "Ifpack_ILUT::Compute: GlobalIndices type unknown";
554 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
558 Xcopy = Teuchos::rcp( &X,
false );
574#ifdef IFPACK_FLOPCOUNTERS
592 const int MaxIters,
const double Tol,
611 if (!
Comm().MyPID()) {
613 os <<
"================================================================================" << endl;
614 os <<
"Ifpack_ILUT: " <<
Label() << endl << endl;
618 os <<
"Relax value = " <<
RelaxValue() << endl;
619 os <<
"Condition number estimate = " <<
Condest() << endl;
625 <<
" % of A)" << endl;
626 os <<
"nonzeros / rows = "
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;
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual int NumProc() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
double ** Pointers() const
int Norm2(double *Result) const
virtual int NumMyRows() const=0
virtual long long NumGlobalNonzeros64() const=0
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual int NumGlobalNonzeros() const=0
virtual long long NumGlobalRows64() 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
bool operator()(const double &x, const double &y)
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 ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
double LevelOfFill_
Level-of-fill.
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
bool IsInitialized_
true if this object has been initialized
double Athresh_
Absolute threshold.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual double InitializeTime() const
Returns the time spent in Initialize().
double DropTolerance() const
Gets the dropping tolerance.
const Epetra_RowMatrix & A_
reference to the matrix to be preconditioned.
int NumApplyInverse_
Contains the number of successful call to 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.
double RelaxValue() const
Set relative threshold value.
virtual int NumCompute() const
Returns the number of calls to Compute().
Teuchos::RefCountPtr< Epetra_Map > SerialMap_
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Teuchos::RefCountPtr< Epetra_SerialComm > SerialComm_
double DropTolerance_
Discards all elements below this tolerance.
double ApplyInverseFlops_
Contain sthe number of flops for 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.
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
U factor.
long long NumGlobalNonzeros64() const
double ComputeTime_
Contains the time for all successful calls to Compute().
virtual int NumInitialize() const
Returns the number of calls to Initialize().
int NumMyRows_
Number of local rows.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
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.
void Destroy()
Releases all allocated memory.
virtual ~Ifpack_ILUT()
Ifpack_ILUT Destructor.
double ComputeFlops_
Contains the number of flops for Compute().
double Relax_
relaxation value
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
double Condest_
Condition number estimate.
double InitializeTime_
Contains the time for all successful calls to Initialize().
int NumCompute_
Contains the number of successful call to Compute().
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
long long GlobalNonzeros_
Global number of nonzeros in L and U factors.
double LevelOfFill() const
Epetra_Time Time_
Used for timing purposed.
double Rthresh_
Relative threshold.
bool IsComputed_
true if this object has been computed
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
L factor.
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
bool UseTranspose_
true if transpose has to be used.
std::string Label_
Label for this object.
double AbsoluteThreshold() const
Get absolute threshold value.
double get(const key_type key)
Returns an element from the hash table, or 0.0 if not found.
void arrayify(key_type *key_array, double *val_array)
Converts the contents in array format for both keys and values.
void set(const key_type key, const double value, const bool addToValue=false)
Sets an element in the hash table.
int getRecommendedHashSize(int n)
void reset()
Resets the entries of the already allocated memory. This method can be used to clean an array,...
int getNumEntries() const
Returns the number of stored entries.