IFPACK Development

Ifpack_CrsIct: A class for constructing and using an incomplete Cholesky factorization of a given Epetra_CrsMatrix. More...
#include <Ifpack_CrsIct.h>
Public Member Functions  
Ifpack_CrsIct (const Epetra_CrsMatrix &A, double Droptol=1.0E4, int Lfil=20)  
Ifpack_CrsIct constuctor with variable number of indices per row.  
Ifpack_CrsIct (const Ifpack_CrsIct &IctOperator)  
Copy constructor.  
virtual  ~Ifpack_CrsIct () 
Ifpack_CrsIct Destructor.  
void  SetAbsoluteThreshold (double Athresh) 
Set absolute threshold value.  
void  SetRelativeThreshold (double Rthresh) 
Set relative threshold value.  
void  SetOverlapMode (Epetra_CombineMode OverlapMode) 
Set overlap mode type.  
int  SetParameters (const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false) 
Set parameters using a Teuchos::ParameterList object.  
int  InitValues (const Epetra_CrsMatrix &A) 
Initialize L and U with values from user matrix A.  
bool  ValuesInitialized () const 
If values have been initialized, this query returns true, otherwise it returns false.  
int  Factor () 
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parameters.  
bool  Factored () const 
If factor is completed, this query returns true, otherwise it returns false.  
int  Solve (bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const 
Returns the result of a Ifpack_CrsIct forward/back solve on a Epetra_MultiVector X in Y.  
int  Multiply (bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const 
Returns the result of multiplying U, D and U^T in that order on an Epetra_MultiVector X in Y.  
int  Condest (bool Trans, double &ConditionNumberEstimate) const 
Returns the maximum over all the condition number estimate for each local ILU set of factors.  
double  GetAbsoluteThreshold () 
Get absolute threshold value.  
double  GetRelativeThreshold () 
Get relative threshold value.  
Epetra_CombineMode  GetOverlapMode () 
Get overlap mode type.  
int  NumGlobalNonzeros () const 
Returns the number of nonzero entries in the global graph.  
int  NumMyNonzeros () const 
Returns the number of nonzero entries in the local graph.  
const Epetra_Vector &  D () const 
Returns the address of the D factor associated with this factored matrix.  
const Epetra_CrsMatrix &  U () const 
Returns the address of the U factor associated with this factored matrix.  
const char *  Label () const 
Returns a character string describing the operator.  
int  SetUseTranspose (bool UseTranspose_in) 
If set true, transpose of this operator will be applied.  
int  Apply (const Epetra_MultiVector &X, Epetra_MultiVector &Y) const 
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.  
int  ApplyInverse (const Epetra_MultiVector &X, Epetra_MultiVector &Y) const 
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.  
double  NormInf () const 
Returns 0.0 because this class cannot compute Infnorm.  
bool  HasNormInf () const 
Returns false because this class cannot compute an Infnorm.  
bool  UseTranspose () const 
Returns the current UseTranspose setting.  
const Epetra_Map &  OperatorDomainMap () const 
Returns the Epetra_Map object associated with the domain of this operator.  
const Epetra_Map &  OperatorRangeMap () const 
Returns the Epetra_Map object associated with the range of this operator.  
const Epetra_Comm &  Comm () const 
Returns the Epetra_BlockMap object associated with the range of this matrix operator.  
Protected Member Functions  
void  SetFactored (bool Flag) 
void  SetValuesInitialized (bool Flag) 
bool  Allocated () const 
int  SetAllocated (bool Flag) 
Friends  
ostream &  operator<< (ostream &os, const Ifpack_CrsIct &A) 
<< operator will work for Ifpack_CrsIct. 
Ifpack_CrsIct: A class for constructing and using an incomplete Cholesky factorization of a given Epetra_CrsMatrix.
The Ifpack_CrsIct class computes a threshold based incomplete LDL^T factorization of a given Epetra_CrsMatrix. The factorization that is produced is a function of several parameters:
Maximum number of entries per row/column in factor  The factorization will contain at most this number of nonzero terms in each row/column of the factorization.
Estimating Preconditioner Condition Numbers
For illconditioned matrices, we often have difficulty computing usable incomplete factorizations. The most common source of problems is that the factorization may encounter a small or zero pivot, in which case the factorization can fail, or even if the factorization succeeds, the factors may be so poorly conditioned that use of them in the iterative phase produces meaningless results. Before we can fix this problem, we must be able to detect it. To this end, we use a simple but effective condition number estimate for .
The condition of a matrix , called , is defined as in some appropriate norm . gives some indication of how many accurate floating point digits can be expected from operations involving the matrix and its inverse. A condition number approaching the accuracy of a given floating point number system, about 15 decimal digits in IEEE double precision, means that any results involving or may be meaningless.
The norm of a vector is defined as the maximum of the absolute values of the vector entries, and the norm of a matrix C is defined as . A crude lower bound for the is where . It is a lower bound because .
For our purposes, we want to estimate , where and are our incomplete factors. Edmond in his Ph.D. thesis demonstrates that provides an effective estimate for . Furthermore, since finding such that is a basic kernel for applying the preconditioner, computing this estimate of is performed by setting , calling the solve kernel to compute and then computing .
A priori Diagonal Perturbations
Given the above method to estimate the conditioning of the incomplete factors, if we detect that our factorization is too illconditioned we can improve the conditioning by perturbing the matrix diagonal and restarting the factorization using this more diagonally dominant matrix. In order to apply perturbation, prior to starting the factorization, we compute a diagonal perturbation of our matrix and perform the factorization on this perturbed matrix. The overhead cost of perturbing the diagonal is minimal since the first step in computing the incomplete factors is to copy the matrix into the memory space for the incomplete factors. We simply compute the perturbed diagonal at this point.
The actual perturbation values we use are the diagonal values with , , where is the matrix dimension and returns the sign of the diagonal entry. This has the effect of forcing the diagonal values to have minimal magnitude of and to increase each by an amount proportional to , and still keep the sign of the original diagonal entry.
Constructing Ifpack_CrsIct objects
Constructing Ifpack_CrsIct objects is a multistep process. The basic steps are as follows:
Note that, even after a matrix is constructed, it is possible to update existing matrix entries. It is not possible to create new entries.
Counting Floating Point Operations
Each Ifpack_CrsIct object keep track of the number of serial floating point operations performed using the specified object as the this argument to the function. The Flops() function returns this number as a double precision number. Using this information, in conjunction with the Epetra_Time class, one can get accurate parallel performance numbers. The ResetFlops() function resets the floating point counter.
Definition at line 159 of file Ifpack_CrsIct.h.
Ifpack_CrsIct::Ifpack_CrsIct  (  const Epetra_CrsMatrix &  A, 
double  Droptol = 1.0E4 , 

int  Lfil = 20 

) 
Ifpack_CrsIct constuctor with variable number of indices per row.
Creates a Ifpack_CrsIct object and allocates storage.
In  A  User matrix to be factored. 
In  Graph  Graph generated by Ifpack_IlukGraph. 
Definition at line 44 of file Ifpack_CrsIct.cpp.
References Zero.
int Ifpack_CrsIct::Apply  (  const Epetra_MultiVector &  X, 
Epetra_MultiVector &  Y  
)  const [inline, virtual] 
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
Note that this implementation of Apply does NOT perform a forward back solve with the LDU factorization. Instead it applies these operators via multiplication with U, D and L respectively. The ApplyInverse() method performs a solve.
In  X  A Epetra_MultiVector of dimension NumVectors to multiply with matrix. 
Out  Y A Epetra_MultiVector of dimension NumVectors containing result. 
Implements Epetra_Operator.
Definition at line 319 of file Ifpack_CrsIct.h.
References UseTranspose().
int Ifpack_CrsIct::ApplyInverse  (  const Epetra_MultiVector &  X, 
Epetra_MultiVector &  Y  
)  const [inline, virtual] 
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
In this implementation, we use several existing attributes to determine how virtual method ApplyInverse() should call the concrete method Solve(). We pass in the UpperTriangular(), the Epetra_CrsMatrix::UseTranspose(), and NoDiagonal() methods. The most notable warning is that if a matrix has no diagonal values we assume that there is an implicit unit diagonal that should be accounted for when doing a triangular solve.
In  X  A Epetra_MultiVector of dimension NumVectors to solve for. 
Out  Y A Epetra_MultiVector of dimension NumVectors containing result. 
Implements Epetra_Operator.
Definition at line 336 of file Ifpack_CrsIct.h.
References UseTranspose().
int Ifpack_CrsIct::Condest  (  bool  Trans, 
double &  ConditionNumberEstimate  
)  const 
Returns the maximum over all the condition number estimate for each local ILU set of factors.
This functions computes a local condition number estimate on each processor and return the maximum over all processor of the estimate.
In  Trans If true, solve transpose problem. 
Out  ConditionNumberEstimate  The maximum across all processors of the infinitynorm estimate of the condition number of the inverse of LDU. 
Definition at line 361 of file Ifpack_CrsIct.cpp.
References Epetra_CrsMatrix::RowMap(), and Solve().
int Ifpack_CrsIct::Factor  (  void  ) 
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parameters.
This function computes the RILU(k) factors L and U using the current:
InitValues() must be called before the factorization can proceed.
Definition at line 227 of file Ifpack_CrsIct.cpp.
References Epetra_CrsMatrix::Comm(), Factored(), Epetra_CrsMatrix::OperatorDomainMap(), Epetra_CrsMatrix::OperatorRangeMap(), Epetra_CrsMatrix::RowMatrixRowMap(), Epetra_Comm::SumAll(), Epetra_CompObject::UpdateFlops(), and View.
int Ifpack_CrsIct::InitValues  (  const Epetra_CrsMatrix &  A  ) 
Initialize L and U with values from user matrix A.
Copies values from the user's matrix into the nonzero pattern of L and U.
In  A  User matrix to be factored. 
Definition at line 150 of file Ifpack_CrsIct.cpp.
References Epetra_CrsMatrix::Comm(), Epetra_Comm::MaxAll(), Epetra_CrsMatrix::OperatorDomainMap(), and Epetra_CrsMatrix::OperatorRangeMap().
int Ifpack_CrsIct::Multiply  (  bool  Trans, 
const Epetra_MultiVector &  X,  
Epetra_MultiVector &  Y  
)  const 
Returns the result of multiplying U, D and U^T in that order on an Epetra_MultiVector X in Y.
In  Trans If true, multiply by L^T, D and U^T in that order. 
In  X  A Epetra_MultiVector of dimension NumVectors to solve for. 
Out  Y A Epetra_MultiVector of dimension NumVectorscontaining result. 
Definition at line 337 of file Ifpack_CrsIct.cpp.
References Epetra_MultiVector::NumVectors(), Epetra_MultiVector::ReciprocalMultiply(), and Epetra_MultiVector::Update().
int Ifpack_CrsIct::SetUseTranspose  (  bool  UseTranspose_in  )  [inline, virtual] 
If set true, transpose of this operator will be applied.
This flag allows the transpose of the given operator to be used implicitly. Setting this flag affects only the Apply() and ApplyInverse() methods. If the implementation of this interface does not support transpose use, this method should return a value of 1.
In  UseTranspose_in If true, multiply by the transpose of operator, otherwise just use operator. 
Implements Epetra_Operator.
Definition at line 305 of file Ifpack_CrsIct.h.
int Ifpack_CrsIct::Solve  (  bool  Trans, 
const Epetra_MultiVector &  X,  
Epetra_MultiVector &  Y  
)  const 
Returns the result of a Ifpack_CrsIct forward/back solve on a Epetra_MultiVector X in Y.
In  Trans If true, solve transpose problem. 
In  X  A Epetra_MultiVector of dimension NumVectors to solve for. 
Out  Y A Epetra_MultiVector of dimension NumVectorscontaining result. 
Definition at line 316 of file Ifpack_CrsIct.cpp.
References Epetra_MultiVector::Multiply(), and Epetra_MultiVector::NumVectors().
Referenced by Condest().