IFPACK Development

IFPACK, developed by Marzio Sala (ETHZ/DINFK) and Micheal Heroux (SNL 9214), provides a suite of objectoriented algebraic preconditioners for the solution of preconditioned iterative solvers. IFPACK constructors expect an Epetra_RowMatrix object for construction. IFPACK is part of the Trilinos Solver Project and IFPACK object interacts well with other Trilinos classes. In particular, IFPACK can be used as a preconditioner for AztecOO.
Most IFPACK preconditioners can be rewritten as additive Schwarz methods, of overlapping domain decomposition type. The user can adopt a minimaloverlap (that is, zerorow overlap), or ask IFPACK to extend the overlap. The resulting preconditioner reads:
\[ P_{IFPACK}^{1} = \sum_{i = 0}^{NumProcs1} P_i A_i^{1} R_i \]
where \(R_i\) is the restriction operator from the global vector, to the overlapping subdomain \(i\), and \(P_i\) is the prolongator operator. \(P_i\) is generally the transpose of \(R_i\) (in which case the resulting preconditioner is symmetic). It is assumed that each subdomain is assigned to a different processor.
A key component of the previous formula is a strategy to apply \(A_i^{1}\). Using IFPACK, this can be defined by one of the following:
More precisely, the choices for \(A_i^{1}\) are:
Point relaxation preconditioners (through class Ifpack_PointRelaxation). Point relaxation preconditioners are probably the simplest iterative methods, and are generally used as smoothers in multilevel methods (for instance, within ML). Available choices are:
Block relaxation preconditioners (through class Ifpack_BlockRelaxation). Block relaxation preconditioners represent an extension of point preconditioners. The local matrix \(A_i\) is divided into blocks, then a relaxation method is applied on the block structure of \(A_i\). Available choices are:
The code extracts the diagonal blocks (the application of whose inverses is required), and stored them either as dense or as sparse matrices. Dense matrices should be used if the dimension of each block is small, since LAPACK routines are used to factorize the block and to perform the backsolve. For large blocks, instead, it is suggested to store the blocks are sparse matrices (Epetra_CrsMatrix's). Any IFPACK preconditioner derived from Ifpack_Preconditioner can be used to apply the inverse of the diagonal blocks.
Point incomplete factorizations (through classes Ifpack_IC, Ifpack_ICT, Ifpack_ILU, Ifpack_ILUT). See General description of incomplete factorizations .
Exact factorizations (through class Ifpack_Amesos):
Chebyshev polynomials (through class Ifpack_Chebyshev).
IFPACK can be downloaded from the web page http://trilinos.sandia.gov/download/
IFPACK is configured with autotools. IFPACK is enabled by default if configured at the Trilinos level (please refer to the Trilinos documentation for more details). This section briefly recalls some of the configure parameters that affect the IFPACK compilation. For a complete list of parameters, please type
% $TRILINOS_HOME/packages/ifpack/configure help
An (incomplete) list of parameters recognized by the IFPACK configure is reported below.
"enableamesos"
enables the support for Amesos (off by default)."enableaztecoo"
enables the support for AztecOO (on by default)."enableteuchos"
enables the support for Teuchos (off by default)."enablethyra"
enables the support for Thyra (off by default)."enabletriutils"
enables the support for Triutils (on by default)."enableifpackmetis"
enables the support for the METIS package. The location of the header files should be specified using "withincdirs=I/include/location"
, the path of the METIS library with "withldflags=L/lib/location"
, and the library name using "withlibs=lmetisname."
.This section details how to use the IFPACK Factory class.
Probably, the easiest way to use IFPACK is through the function class Ifpack. This is a factory class, that contains only one method, Create(). A call to Create() returns a pointer to a newlycreated IFPACK object. The user is responsible of deleting this object.
An example of usage is reported by the following code:
#include "Ifpack.h" ... Ifpack Factory; Epetra_RowMatrix* A; // A is already FillComplete()'d // but its values can still change (not its structure) string PrecType = "Amesos"; // exact solve on each subdomain int OverlapLevel = 1; // one row of overlap among the processes Ifpack_Preconditioner* Prec = Factory.Create(PrecType, A, OverlapLevel); assert (Prec != 0); IFPACK_CHK_ERR(Prec>SetParameters(List)); IFPACK_CHK_ERR(Prec>Initialize()); IFPACK_CHK_ERR(Prec>Compute()); // we can query the status of the preconditioner assert (Prec>IsInitialized() == true); assert (Prec>IsComputed() == true); // information about the preconditioner // 1. number of initialization phases cout << Prec>NumInitialize() << endl; // 2. number of computation phases cout << Prec>NumCompute() << endl; // 3. number of applications of the prconditioner cout << Prec>NumApplyInverse() << endl; // we can query for flops and CPU time as well, consult the // documentation of Ifpack_Preconditioner, or simply write cout << *Prec; // now destroy the preconditioner delete Prec;
Although the Ifpack factory is appropriate for most users, it is possible to create IFPACK preconditioners by handling directly an Ifpack_AdditiveSchwarz object. The preconditioner of the previous example is equivalent to the one created in the following fragment of code:
#include "Ifpack_AdditiveSchwarz.h" ... Epetra_RowMatrix* A; int OverlapLevel = 1; Ifpack_AdditiveSchwarz Prec<Ifpack_Amesos> (A, OverlapLevel);
Another example is the following, where the local solve is a block Jacobi method, which uses dense containers (LAPACK).
#include "Ifpack_AdditiveSchwarz.h" #include "Ifpack_DenseContainer.h" #include "Ifpack_BlockJacobi.h" ... Epetra_RowMatrix* A; int OverlapLevel = 1; Ifpack_AdditiveSchwarz Prec<Ifpack_BlockJacobi<Ifpack_DenseContainer> > (A, OverlapLevel);
Using sparse containers and Amesos the code may look as follows.
#include "Ifpack_AdditiveSchwarz.h" #include "Ifpack_SparseContainer.h" #include "Ifpack_BlockJacobi.h" #include "Ifpack_Amesos.h" ... Epetra_RowMatrix* A; int OverlapLevel = 1; Ifpack_AdditiveSchwarz Prec<Ifpack_BlockJacobi<Ifpack_SparseContainer<Ifpack_Amesos> > > (A, OverlapLevel);
Other examples are reported in the example
subdirectory with several comments:
The complete list of supported parameters is reported below.
"relaxation: type"
: specifies the type of point and block relaxation scheme. Valid values:"Jacobi"
;"GaussSeidel"
;"symmetric GaussSeidel"
;"relaxation: sweeps"
: specifies the number of sweeps in the application of point relaxation schemes."relaxation: damping factor"
: specifies the damping factor in the application of point relaxation schemes."relaxation: min diagonal value"
: replace diagonal values below this value with this value (only for point relaxation)."relaxation: zero starting solution"
: if true
, uses the values in the vector to be preconditioned as starting solution. Otherwise, use the zero vector as starting solution."relaxation: backward mode"
: Does GaussSeidel in backward order (rather than forward)"relaxation: use l1"
: Use the "L1GaussSeidel" method of Baker, Falgout, Kolev and Yang (SISC 2011). This can be especially advantageous for using GaussSeidel and symmetric GaussSeidel in parallel."relaxation: l1 eta"
: Sets the "eta" parameter used by the L1 GS/SGS method."relaxation: number of local smoothing indices"
: Sets the number of local unknowns over which to smooth. This is used for selective (local) smoothing."relaxation: local smoothing indices"
: A list of local unknowns over which to smooth. This can be used for selective (local) smoothing or to force the method to smooth unknowns in a particular order (such as upwind smoothing)."partitioner: type"
: specifies the scheme to adopt to partition the graph (for block relaxation methods only). Valid values:"linear"
(uses a linear decomposition);"greedy"
(uses a simple greedy algorithm);"metis"
(calls METIS, this required enableifpackmetis
)."equation"
(groups together all the nodes belonging to the same equation. The number of equation is specified using "partitioner: local
parts"
, and it is supposed that the equation for each grid node are ordered consecutevely)."partitioner: overlap"
: specifies the overlap among the blocks which can differ from the overlap among the processors. Note that only the Jacobi block relaxation scheme can take advantage of nonzero overlaps."partitioner: local parts"
: specifies the number of local blocks."partitioner: root node"
: specifies the root node for greedy algorithm."partitioner: use symmetric graph"
: if true
, METIS will partition the graph of A + A^{T}. If false
, METIS will partition the graph of A. Note that METIS can core dump if the input graph is nonsymmetric. Users should set this option to false
only when the graph is symmetric. If the nonsymmetry is determined by Dirichlet nodes, then the singleton filter should create a symmetric graph. Note also that dropping techniques applied to nonsymmetric matrices can result in nonsymmetric graph."amesos: solver type"
: defines the Amesos solver to be used by class Ifpack_Amesos. Valid values:"Amesos_Lapack"
(use Amesos interface to LAPACK);"Amesos_Klu"
(use Amesos' internal solver KLU);"Amesos_Umfpack"
(use UMFPACK interface);"Amesos_Superlu"
(use serial SuperLU interface);"Amesos_Mumps"
(use MUMPS interface);"Amesos_Dscpack"
(use DSCPACK interface)."fact: leveloffill"
: defines the level of fill for ILU factorizations. This is based on powers of the graph, so the value 0 means nofill."fact: ilut leveloffill"
: defines the level of fill for ILUT factorizations. This is a ratio, so 1.0 means same number of nonzeros for the ILU factors as in the original matrix."fact: ict leveloffill"
: defines the level of fill for ICT factorizations, and is a fillratio as for ILUT. Currently this parameter is used in both Ifpack_IC and Ifpack_ICT."fact: absolute threshold"
: defines the value \(\alpha\) to add to each diagonal element (times the signum of the actual diagonal element), for incomplete factorizations only."fact: relative threshold"
: defines \(\rho\), so that the diagonal element of the matrix (without the contribution specified by "fact: absolute threshold"
) is multiplied by \(rho \)."fact: relax value"
: if different from zero, the elements dropped during the factorization process will be added to the diagonal term, multiplied by the specified value."schwarz: combine mode"
: defines how values corresponding to overlapping nodes are handled (after the solution of the local problem). Any Epetra_CombineMode is valid. Users should set this value to Add
if interested in a symmetric preconditioner. Otherwise, the default value of "Zero"
usually results in better convergence. Other details on combine mode are reported. All Ifpack preconditioners work on parallel distributed memory computers by using the row partitioning the user input matrix to determine the partitioning for local ILU 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 selfcontained local matrix and all column entries that reach to offprocessor 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 ILU factorization process. Level of overlap is defined during the construction of the Ifpack_IlukGraph object. The application of the inverse of the local matrix results in redundant approximations for any elements of y that correspond to rows that are part of more than one local ILU factor. The combine mode defines how these redundancies are handled using the Epetra_CombineMode enum. The default (Zero
) means to zero out all values of y for rows that were not part of the original matrix row distribution."schwarz: compute condest"
: it true
, compute a cheap estimate of the condition number. Warning: this is not the REAL conditioner number, but only a very lightweight indication of how well the current preconditioner is expected to behave."schwarz: filter singletons"
: it true
, singletons in the local matrix are removed."schwarz: reordering type"
: it none
, no reordering is applied. If rcm
or metis
, a reverse CuthillMcKee or METIS are used to reorder to local matrix.The table below reports the C++ type and the default value of each parameter.
Parameter name  C++ type  Defaut value 
"relaxation: type"  string  "Jacobi" 
"relaxation: sweeps"  int  "1" 
"relaxation: damping factor"  double  "1.0" 
"relaxation: min diagonal value"  double  "0.0" 
"relaxation: zero starting solution"  bool  "true" 
"relaxation: backward mode"  bool  "false" 
"relaxation: use l1"  bool  "false" 
"relaxation: l1 eta"  double  "1.5" 
"relaxation: number of local smoothing indices"  int  "A>NumMyRows()" 
"relaxation: local smoothing indices"  int*  "0" 
"partitioner: type"  string  "greedy" 
"partitioner: overlap"  int  "0" 
"partitioner: local parts"  int  "1" 
"partitioner: root node"  int  "0" 
"partitioner: use symmetric graph"  bool  "true" 
"amesos: solver type"  string  "Amesos_Klu" 
"fact: absolute threshold"  double  "0.0" 
"fact: relative threshold"  double  "0.0" 
"fact: relax value"  double  "0.0" 
"schwarz: combine mode"  Epetra_CombineMode  "Zero" 
"schwarz: compute condest"  bool  "true" 
"schwarz: reordering type"  string  "none" 
"schwarz: filter singletons"  string  "false" 
IFPACK classes are divided into classes for users and classes for developers.
Users' Tools
Developers' Tools
The following table defines the class of IFPACK errors. Return code should be checked using the IFPACK_CHK_ERR() macro. In general terms, we follow this convention:
Error Code  Meaning 
1  Generic Error (called method or function returned an error) 
2  Input data not valid (wrong parameter, outofbounds, wrong dimensions, matrix is not square,...) 
3  Data has not been correctly preprocessed 
4  Problem encountered during application of the algorithm (division by zero, outofbounds, ...) 
5  Memory allocation error 
98  Feature is not supported 
99  Feature is not implemented yet (check Known Bugs and Future Developments, or submit a bug) 
UseTranspose()
==
true. Ifpack_PointRelaxation and Ifpack_BlockRelaxation do not support UseTranspose
, while Ifpack_ILU, Ifpack_ILUT and Ifpack_Amesos do.Under terms of Contract DEAC0494AL85000, there is a nonexclusive license for use of this work by or on behalf of the U.S. Government. This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 021101301 USA
You can browse all of Ifpack as a single doxygen collection. Warning, this is not the recommended way to learn about Ifpack software. However, this is a good way to browse the directory structure of ifpack, to locate files, etc.