The PCG (Preconditioned Conjugate Gradient [HS1952] implementation of
the SUNLinearSolver
module provided with SUNDIALS, SUNLinSol_PCG,
is an iterative linear solver that is designed to be compatible with
any N_Vector
implementation (serial, threaded, parallel, and
user-supplied) that supports a minimal subset of operations
(N_VClone()
, N_VDotProd()
, N_VScale()
,
N_VLinearSum()
, N_VProd()
, and
N_VDestroy()
). Unlike the SPGMR and SPFGMR algorithms, PCG
requires a fixed amount of memory that does not increase with the
number of allowed iterations.
Unlike all of the other iterative linear solvers supplied with SUNDIALS, PCG should only be used on symmetric linear systems (e.g. mass matrix linear systems encountered in ARKode). As a result, the explanation of the role of scaling and preconditioning matrices given in general must be modified in this scenario. The PCG algorithm solves a linear system \(Ax = b\) where \(A\) is a symmetric (\(A^T=A\)), real-valued matrix. Preconditioning is allowed, and is applied in a symmetric fashion on both the right and left. Scaling is also allowed and is applied symmetrically. We denote the preconditioner and scaling matrices as follows:
The matrices \(A\) and \(P\) are not required explicitly; only routines
that provide \(A\) and \(P^{-1}\) as operators are required. The diagonal
of the matrix \(S\) is held in a single N_Vector
, supplied by the user.
In this notation, PCG applies the underlying CG algorithm to the equivalent transformed system
where
The scaling matrix must be chosen so that the vectors \(SP^{-1}b\) and \(S^{-1}Px\) have dimensionless components.
The stopping test for the PCG iterations is on the L2 norm of the scaled preconditioned residual:
where \(\| v \|_S = \sqrt{v^T S^T S v}\), with an input tolerance \(\delta\).
The header file to be included when using this module
is sunlinsol/sunlinsol_pcg.h
. The SUNLinSol_PCG module
is accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsolpcg
module library.
The module SUNLinSol_PCG provides the following user-callable routines:
SUNLinSol_PCG
(N_Vector y, int pretype, int maxl)¶This constructor function creates and allocates memory for a PCG
SUNLinearSolver
. Its arguments are an N_Vector
, a flag
indicating to use preconditioning, and the number of linear
iterations to allow.
This routine will perform consistency checks to ensure that it is
called with a consistent N_Vector
implementation (i.e. that it
supplies the requisite vector operations). If y
is
incompatible then this routine will return NULL
.
A maxl
argument that is \(\le0\) will result in the default
value (5).
Since the PCG algorithm is designed to only support symmetric
preconditioning, then any of the pretype
inputs PREC_LEFT
(1), PREC_RIGHT
(2), or PREC_BOTH
(3) will result in use
of the symmetric preconditioner; any other integer input will
result in the default (no preconditioning). Although some SUNDIALS
solvers are designed to only work with left preconditioning (IDA
and IDAS) and others with only right preconditioning (KINSOL), PCG
should only be used with these packages when the linear systems
are known to be symmetric. Since the scaling of matrix rows and
columns must be identical in a symmetric matrix, symmetric
preconditioning should work appropriately even for packages
designed with one-sided preconditioning in mind.
SUNLinSol_PCGSetPrecType
(SUNLinearSolver S, int pretype)¶This function updates the flag indicating use of preconditioning.
As above, any one of the input values, PREC_LEFT
(1),
PREC_RIGHT
(2), or PREC_BOTH
(3) will enable
preconditioning; PREC_NONE
(0) disables preconditioning.
This routine will return with one of the error codes
SUNLS_ILL_INPUT
(illegal pretype
), SUNLS_MEM_NULL
(S
is NULL
), or SUNLS_SUCCESS
.
SUNLinSol_PCGSetMaxl
(SUNLinearSolver S, int maxl)¶This function updates the number of linear solver iterations to allow.
A maxl
argument that is \(\le0\) will result in the default
value (5).
This routine will return with one of the error codes
SUNLS_MEM_NULL
(S
is NULL
) or SUNLS_SUCCESS
.
SUNLinSolSetInfoFile_PCG
(SUNLinearSolver LS, FILE* info_file)¶The function SUNLinSolSetInfoFile_PCG()
sets the
output file where all informative (non-error) messages should be directed.
stdout
by default);NULL
input will disable outputNULL
Notes:
This function is intended for users that wish to monitor the linear
solver progress. By default, the file pointer is set to stdout
.
SUNDIALS must be built with the CMake option ``SUNDIALS_BUILD_WITH_MONITORING``, to utilize this function. See section Configuration options (Unix/Linux) for more information.
SUNLinSolSetPrintLevel_PCG
(SUNLinearSolver LS, int print_level)¶The function SUNLinSolSetPrintLevel_PCG()
specifies the
level of verbosity of the output.
LS – a SUNLinSol object
print_level – flag indicating level of verbosity; must be one of:
- 0, no information is printed (default)
- 1, for each linear iteration the residual norm is printed
NULL
Notes: This function is intended for users that wish to monitor the linear solver progress. By default, the print level is 0.
SUNDIALS must be built with the CMake option ``SUNDIALS_BUILD_WITH_MONITORING``, to utilize this function. See section Configuration options (Unix/Linux) for more information.
For backwards compatibility, we also provide the wrapper functions, each with identical input and output arguments to the routines that they wrap:
SUNPCG
(N_Vector y, int pretype, int maxl)¶Wrapper function for SUNLinSol_PCG()
SUNPCGSetPrecType
(SUNLinearSolver S, int pretype)¶Wrapper function for SUNLinSol_PCGSetPrecType()
SUNPCGSetMaxl
(SUNLinearSolver S, int maxl)¶Wrapper function for SUNLinSol_PCGSetMaxl()
For solvers that include a Fortran interface module, the
SUNLinSol_PCG module also includes the Fortran-callable
function FSUNPCGInit()
to initialize
this SUNLinSol_PCG module for a given SUNDIALS solver.
FSUNPCGInit
(CODE, PRETYPE, MAXL, IER)¶Initializes a PCG SUNLinearSolver
structure for
use in a SUNDIALS package.
This routine must be called after the N_Vector
object has
been initialized.
int
, input) – flag denoting the SUNDIALS solver
this matrix will be used for: CVODE=1, IDA=2, KINSOL=3, ARKode=4.int
, input) – flag denoting whether to use
symmetric preconditioning: no=0, yes=1.int
, input) – number of PCG iterations to allow.int
, output) – return flag (0 success, -1 for failure).Additionally, when using ARKode with a non-identity
mass matrix, the Fortran-callable function
FSUNMassPCGInit()
initializes this
SUNLinSol_PCG module for solving mass matrix linear systems.
FSUNMassPCGInit
(PRETYPE, MAXL, IER)¶Initializes a PCG SUNLinearSolver
structure for
use in solving mass matrix systems in ARKode.
This routine must be called after the N_Vector
object has
been initialized.
int
, input) – flag denoting whether to use
symmetric preconditioning: no=0, yes=1.int
, input) – number of PCG iterations to allow.int
, output) – return flag (0 success, -1 for failure).The SUNLinSol_PCGSetPrecType()
and SUNLinSol_PCGSetMaxl()
routines also support Fortran interfaces for the system and mass
matrix solvers:
FSUNPCGSetPrecType
(CODE, PRETYPE, IER)¶Fortran interface to SUNLinSol_PCGSetPrecType()
for system
linear solvers.
This routine must be called after FSUNPCGInit()
has
been called.
Arguments: all should have type int
, and have meanings
identical to those listed above.
FSUNMassPCGSetPrecType
(PRETYPE, IER)¶Fortran interface to SUNLinSol_PCGSetPrecType()
for mass matrix
linear solvers in ARKode.
This routine must be called after FSUNMassPCGInit()
has
been called.
Arguments: all should have type int
, and have meanings
identical to those listed above.
FSUNPCGSetMaxl
(CODE, MAXL, IER)¶Fortran interface to SUNLinSol_PCGSetMaxl()
for system
linear solvers.
This routine must be called after FSUNPCGInit()
has
been called.
Arguments: all should have type int
, and have meanings
identical to those listed above.
FSUNMassPCGSetMaxl
(MAXL, IER)¶Fortran interface to SUNLinSol_PCGSetMaxl()
for mass matrix
linear solvers in ARKode.
This routine must be called after FSUNMassPCGInit()
has
been called.
Arguments: all should have type int
, and have meanings
identical to those listed above.
The SUNLinSol_PCG module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_PCG {
int maxl;
int pretype;
int numiters;
realtype resnorm;
int last_flag;
ATimesFn ATimes;
void* ATData;
PSetupFn Psetup;
PSolveFn Psolve;
void* PData;
N_Vector s;
N_Vector r;
N_Vector p;
N_Vector z;
N_Vector Ap;
int print_level;
FILE* info_file;
};
These entries of the content field contain the following information:
maxl
- number of PCG iterations to allow (default is 5),pretype
- flag for use of preconditioning (default is none),numiters
- number of iterations from the most-recent solve,resnorm
- final linear residual norm from the most-recent
solve,last_flag
- last error return flag from an internal
function,ATimes
- function pointer to perform \(Av\) product,ATData
- pointer to structure for ATimes
,Psetup
- function pointer to preconditioner setup routine,Psolve
- function pointer to preconditioner solve routine,PData
- pointer to structure for Psetup
and Psolve
,s
- vector pointer for supplied scaling matrix
(default is NULL
),r
- a N_Vector
which holds the preconditioned linear system
residual,p, z, Ap
- N_Vector
used for workspace by the
PCG algorithm.print_level
- controls the amount of information to be printed to the info fileinfo_file
- the file where all informative (non-error) messages will be directedThis solver is constructed to perform the following operations:
N_Vector
solver data is allocated, with
vectors cloned from a template N_Vector
that is input, and
default solver parameters are set.ATimes
, PSetup
, and Psolve
function pointers and
s
scaling vector.NULL
PSetup
function is
called. Typically, this is provided by the SUNDIALS solver
itself, that translates between the generic PSetup
function and
the solver-specific routine (solver-supplied or user-supplied).The SUNLinSol_PCG module defines implementations of all “iterative” linear solver operations listed in the section The SUNLinearSolver API:
SUNLinSolGetType_PCG
SUNLinSolInitialize_PCG
SUNLinSolSetATimes_PCG
SUNLinSolSetPreconditioner_PCG
SUNLinSolSetScalingVectors_PCG
– since PCG only supports
symmetric scaling, the second N_Vector
argument to this function
is ignoredSUNLinSolSetup_PCG
SUNLinSolSolve_PCG
SUNLinSolNumIters_PCG
SUNLinSolResNorm_PCG
SUNLinSolResid_PCG
SUNLinSolLastFlag_PCG
SUNLinSolSpace_PCG
SUNLinSolFree_PCG