IFPACK Development
Public Member Functions | Protected Member Functions | Friends
Ifpack_CrsRick Class Reference

Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a given Epetra_CrsMatrix. More...

#include <Ifpack_CrsRick.h>

Inheritance diagram for Ifpack_CrsRick:
Inheritance graph
[legend]
Collaboration diagram for Ifpack_CrsRick:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 Ifpack_CrsRick (const Epetra_CrsMatrix &A, const Ifpack_IlukGraph &Graph)
 Ifpack_CrsRick constuctor with variable number of indices per row.
 Ifpack_CrsRick (const Ifpack_CrsRick &Matrix)
 Copy constructor.
virtual ~Ifpack_CrsRick ()
 Ifpack_CrsRick Destructor.
int InitValues ()
 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.
void SetRelaxValue (double RelaxValue)
 Set RIC(k) relaxation parameter.
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 &parameterlist, bool cerr_warning_if_unused=false)
 Set parameters using a Teuchos::ParameterList object.
int Factor ()
 Compute IC factor L 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_Vector &x, Epetra_Vector &y) const
 Returns the result of a Ifpack_CrsRick forward/back solve on a Epetra_Vector x in y.
int Solve (bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
 Returns the result of a Ifpack_CrsRick 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 IC set of factors.
double GetRelaxValue ()
 Get RIC(k) relaxation parameter.
double GetAbsoluteThreshold ()
 Get absolute threshold value.
double GetRelativeThreshold ()
 Get relative threshold value.
Epetra_CombineMode GetOverlapMode ()
 Get overlap mode type.
int NumGlobalRows () const
 Returns the number of global matrix rows.
int NumGlobalCols () const
 Returns the number of global matrix columns.
int NumGlobalNonzeros () const
 Returns the number of nonzero entries in the global graph.
virtual int NumGlobalDiagonals () const
 Returns the number of diagonal entries found in the global input graph.
int NumMyRows () const
 Returns the number of local matrix rows.
int NumMyCols () const
 Returns the number of local matrix columns.
int NumMyNonzeros () const
 Returns the number of nonzero entries in the local graph.
virtual int NumMyDiagonals () const
 Returns the number of diagonal entries found in the local input graph.
int IndexBase () const
 Returns the index base for row and column indices for this graph.
const Ifpack_IlukGraphGraph () const
 Returns the address of the Ifpack_IlukGraph associated with this factored matrix.
const Epetra_VectorD () const
 Returns the address of the D factor associated with this factored matrix.
const Epetra_CrsMatrixU () 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)
 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 Inf-norm.
bool HasNormInf () const
 Returns false because this class cannot compute an Inf-norm.
bool UseTranspose () const
 Returns the current UseTranspose setting.
const Epetra_MapOperatorDomainMap () const
 Returns the Epetra_Map object associated with the domain of this operator.
const Epetra_MapOperatorRangeMap () const
 Returns the Epetra_Map object associated with the range of this 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_CrsRick &A)
 << operator will work for Ifpack_CrsRick.

Detailed Description

Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a given Epetra_CrsMatrix.

The Ifpack_CrsRick class computes a "Relaxed" IC factorization with level k fill of a given Epetra_CrsMatrix. The factorization that is produced is a function of several parameters:

  1. The pattern of the matrix - All fill is derived from the original matrix nonzero structure. Level zero fill is defined as the original matrix pattern (nonzero structure), even if the matrix value at an entry is stored as a zero. (Thus it is possible to add entries to the IC factors by adding zero entries the original matrix.)

  2. Level of fill - Starting with the original matrix pattern as level fill of zero, the next level of fill is determined by analyzing the graph of the previous level and determining nonzero fill that is a result of combining entries that were from previous level only (not the current level). This rule limits fill to entries that are direct decendents from the previous level graph. Fill for level k is determined by applying this rule recursively. For sufficiently large values of k, the fill would eventually be complete and an exact Cholesky factorization would be computed. Level of fill is defined during the construction of the Ifpack_IlukGraph object.

  3. Level of overlap - All Ifpack preconditioners work on parallel distributed memory computers by using the row partitioning the user input matrix to determine the partitioning for local IC factors. If the level of overlap is set to zero, the rows of the user matrix that are stored on a given processor are treated as a self-contained local matrix and all column entries that reach to off-processor entries are ignored. Setting the level of overlap to one tells Ifpack to increase the size of the local matrix by adding rows that are reached to by rows owned by this processor. Increasing levels of overlap are defined recursively in the same way. For sufficiently large levels of overlap, the entire matrix would be part of each processor's local IC factorization process. Level of overlap is defined during the construction of the Ifpack_IlukGraph object.

    Once the factorization is computed, applying the factorization \(LL^{T}y = x\) results in redundant approximations for any elements of y that correspond to rows that are part of more than one local IC factor. The OverlapMode (changed by calling SetOverlapMode()) defines how these redundancies are handled using the Epetra_CombineMode enum. The default is to zero out all values of y for rows that were not part of the original matrix row distribution.

  4. Fraction of relaxation - Ifpack_CrsRick computes the IC factorization row-by-row. As entries at a given row are computed, some number of them will be dropped because they do match the prescribed sparsity pattern. The relaxation factor determines how these dropped values will be handled. If the RelaxValue (changed by calling SetRelaxValue()) is zero, then these extra entries will by dropped. This is a classical IC approach. If the RelaxValue is 1, then the sum of the extra entries will be added to the diagonal. This is a classical Modified IC (MIC) approach. If RelaxValue is between 0 and 1, then RelaxValue times the sum of extra entries will be added to the diagonal.

    For most situations, RelaxValue should be set to zero. For certain kinds of problems, e.g., reservoir modeling, there is a conservation principle involved such that any operator should obey a zero row-sum property. MIC was designed for these cases and you should set the RelaxValue to 1. For other situations, setting RelaxValue to some nonzero value may improve the stability of factorization, and can be used if the computed IC factors are poorly conditioned.

  5. Diagonal perturbation - Prior to computing the factorization, it is possible to modify the diagonal entries of the matrix for which the factorization will be computing. If the absolute and relative perturbation values are zero and one, respectively, the factorization will be compute for the original user matrix A. Otherwise, the factorization will computed for a matrix that differs from the original user matrix in the diagonal values only. Below we discuss the details of diagonal perturbations. The absolute and relative threshold values are set by calling SetAbsoluteThreshold() and SetRelativeThreshold(), respectively.

Estimating Preconditioner Condition Numbers

For ill-conditioned 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 $(LU)^{-1}$.

The condition of a matrix $B$, called $cond_p(B)$, is defined as $cond_p(B) = \|B\|_p\|B^{-1}\|_p$ in some appropriate norm $p$. $cond_p(B)$ 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 $B$ or $B^{-1}$ may be meaningless.

The $\infty$-norm of a vector $y$ is defined as the maximum of the absolute values of the vector entries, and the $\infty$-norm of a matrix C is defined as $\|C\|_\infty = \max_{\|y\|_\infty = 1} \|Cy\|_\infty$. A crude lower bound for the $cond_\infty(C)$ is $\|C^{-1}e\|_\infty$ where $e = (1, 1, \ldots, 1)^T$. It is a lower bound because $cond_\infty(C) = \|C\|_\infty\|C^{-1}\|_\infty \ge \|C^{-1}\|_\infty \ge |C^{-1}e\|_\infty$.

For our purposes, we want to estimate $cond_\infty(LU)$, where $L$ and $U$ are our incomplete factors. Edmond in his Ph.D. thesis demonstrates that $\|(LU)^{-1}e\|_\infty$ provides an effective estimate for $cond_\infty(LU)$. Furthermore, since finding $z$ such that $LUz = y$ is a basic kernel for applying the preconditioner, computing this estimate of $cond_\infty(LU)$ is performed by setting $y = e$, calling the solve kernel to compute $z$ and then computing $\|z\|_\infty$.

A priori Diagonal Perturbations

Given the above method to estimate the conditioning of the incomplete factors, if we detect that our factorization is too ill-conditioned 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 $A$ 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 $A$ 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 $(d_1, d_2, \ldots, d_n)$ with $d_i = sgn(d_i)\alpha + d_i\rho$, $i=1, 2, \ldots, n$, where $n$ is the matrix dimension and $sgn(d_i)$ returns the sign of the diagonal entry. This has the effect of forcing the diagonal values to have minimal magnitude of $\alpha$ and to increase each by an amount proportional to $\rho$, and still keep the sign of the original diagonal entry.

Constructing Ifpack_CrsRick objects

Constructing Ifpack_CrsRick objects is a multi-step process. The basic steps are as follows:

  1. Create Ifpack_CrsRick instance, including storage, via constructor.
  2. Enter values via one or more Put or SumInto functions.
  3. Complete construction via FillComplete call.

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_CrsRick 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.

Warning:
A Epetra_Map is required for the Ifpack_CrsRick constructor.

Definition at line 192 of file Ifpack_CrsRick.h.


Constructor & Destructor Documentation

Ifpack_CrsRick::Ifpack_CrsRick ( const Epetra_CrsMatrix A,
const Ifpack_IlukGraph Graph 
)

Ifpack_CrsRick constuctor with variable number of indices per row.

Creates a Ifpack_CrsRick object and allocates storage.

Parameters:
InA - User matrix to be factored.
InGraph - Graph generated by Ifpack_IlukGraph.

Definition at line 42 of file Ifpack_CrsRick.cpp.

References Zero.


Member Function Documentation

int Ifpack_CrsRick::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.

Parameters:
InX - A Epetra_MultiVector of dimension NumVectors to multiply with matrix.
OutY -A Epetra_MultiVector of dimension NumVectors containing result.
Returns:
Integer error code, set to 0 if successful.

Implements Epetra_Operator.

Definition at line 392 of file Ifpack_CrsRick.h.

References Multiply(), and UseTranspose().

int Ifpack_CrsRick::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.

Parameters:
InX - A Epetra_MultiVector of dimension NumVectors to solve for.
OutY -A Epetra_MultiVector of dimension NumVectors containing result.
Returns:
Integer error code, set to 0 if successful.

Implements Epetra_Operator.

Definition at line 409 of file Ifpack_CrsRick.h.

References Solve(), and UseTranspose().

int Ifpack_CrsRick::Condest ( bool  Trans,
double &  ConditionNumberEstimate 
) const

Returns the maximum over all the condition number estimate for each local IC set of factors.

This functions computes a local condition number estimate on each processor and return the maximum over all processor of the estimate.

Parameters:
InTrans -If true, solve transpose problem.
OutConditionNumberEstimate - The maximum across all processors of the infinity-norm estimate of the condition number of the inverse of LDU.

Definition at line 551 of file Ifpack_CrsRick.cpp.

References Epetra_CrsMatrix::RowMap(), and Solve().

int Ifpack_CrsRick::Factor ( void  )

Compute IC factor L using the specified graph, diagonal perturbation thresholds and relaxation parameters.

This function computes the RIC(k) factor L using the current:

  1. Ifpack_IlukGraph specifying the structure of L.
  2. Value for the RIC(k) relaxation parameter.
  3. Value for the a priori diagonal threshold values.

InitValues() must be called before the factorization can proceed.

Definition at line 221 of file Ifpack_CrsRick.cpp.

References Epetra_BlockMap::Comm(), Epetra_CrsMatrix::ExtractMyRowCopy(), Epetra_CrsMatrix::ExtractMyRowView(), Factored(), Ifpack_IlukGraph::L_Graph(), Epetra_CrsMatrix::MaxNumEntries(), NumMyCols(), NumMyRows(), Epetra_CrsMatrix::ReplaceMyValues(), Epetra_CrsGraph::RowMap(), Epetra_Comm::SumAll(), Epetra_CompObject::UpdateFlops(), and ValuesInitialized().

int Ifpack_CrsRick::InitValues ( )
int Ifpack_CrsRick::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.

Parameters:
InTrans -If true, multiply by L^T, D and U^T in that order.
InX - A Epetra_MultiVector of dimension NumVectors to solve for.
OutY -A Epetra_MultiVector of dimension NumVectorscontaining result.
Returns:
Integer error code, set to 0 if successful.

Definition at line 486 of file Ifpack_CrsRick.cpp.

References Epetra_BlockMap::DistributedGlobal(), Epetra_CrsGraph::DomainMap(), Epetra_DistObject::Export(), Epetra_CompObject::GetFlopCounter(), Epetra_DistObject::Import(), Insert, Ifpack_IlukGraph::L_Graph(), Ifpack_IlukGraph::LevelOverlap(), Epetra_CrsMatrix::Multiply(), Epetra_MultiVector::NumVectors(), Ifpack_IlukGraph::OverlapGraph(), Ifpack_IlukGraph::OverlapImporter(), Epetra_MultiVector::ReciprocalMultiply(), Epetra_CrsGraph::RowMap(), Epetra_CompObject::SetFlopCounter(), and Epetra_MultiVector::Update().

Referenced by Apply().

int Ifpack_CrsRick::SetUseTranspose ( bool  UseTranspose) [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.

Parameters:
InUseTranspose -If true, multiply by the transpose of operator, otherwise just use operator.
Returns:
Always returns 0.

Implements Epetra_Operator.

Definition at line 378 of file Ifpack_CrsRick.h.

References UseTranspose().

int Ifpack_CrsRick::Solve ( bool  Trans,
const Epetra_MultiVector X,
Epetra_MultiVector Y 
) const
int Ifpack_CrsRick::Solve ( bool  Trans,
const Epetra_Vector x,
Epetra_Vector y 
) const

The documentation for this class was generated from the following files:
 All Classes Files Functions Variables Enumerations Friends