Most IFPACK preconditioners can be rewritten as additive Schwarz methods, of overlapping domain decomposition type. The user can adopt a minimal-overlap (that is, zero-row overlap), or ask IFPACK to extend the overlap. The resulting preconditioner reads:
where is the restriction operator from the global vector, to the overlapping subdomain , and is the prolongator operator. is generally be the transpose of (in which case the resulting preconditioner is symmetic). If is supposed that each subdomain is assigned to a different processors.
A key component of the previous formula is a strategy to apply . Using IFPACK, this can be defined by one of the following:
More precisely, the choices for 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.
% $TRILINOS_HOME/packages/ifpack/configure --help
An (incomplete) list of parameters recognized by the IFPACK configure is reported below.
"enable-amesos"enables the support for Amesos (off by default).
"enable-aztecoo"enables the support for AztecOO (on by default).
"enable-teuchos"enables the support for Teuchos (off by default).
"enable-triutils"enables the support for Triutils (on by default).
"enable-ifpack-metis"enables the support for the METIS package. The location of the header files should be specified using
"--with-incdirs=-I/include/location", the path of the METIS library with
"--with-ldflags=-L/lib/location", and the library name using
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 newly-created 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:
"relaxation: type": specifies the type of point and block relaxation scheme. Valid values:
"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.
"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
"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 non-zero 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 + AT. If
false, METIS will partition the graph of A. Note that METIS can core dump if the input graph is non-symmetric. Users should set this option to
falseonly when the graph is symmetric. If the non-symmetry is determined by Dirichlet nodes, then the singleton filter should create a symmetric graph. Note also that dropping techniques applied to non-symmetric matrices can result in non-symmetric 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: level-of-fill": defines the level of fill for IC and ILU factorizations.
"fact: ilut level-of-fill": defines the level of fill for ICT and ILUT factorizations.
"fact: absolute threshold": defines the value to add to each diagonal element (times the signum of the actual diagonal element), for incomplete factorizations only.
"fact: relative threshold": defines , so that the diagonal element of the matrix (without the contribution specified by
"fact: absolute threshold") is multiplied by .
"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
Addif 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 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 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 light-weight 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
metis, a reverse Cuthill-McKee 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|
|-1||Generic Error (called method or function returned an error)|
|-2||Input data not valid (wrong parameter, out-of-bounds, wrong dimensions, matrix is not square,...)|
|-3||Data has not been correctly pre-processed|
|-4||Problem encountered during application of the algorithm (division by zero, out-of-bounds, ...)|
|-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)|
Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA