123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285 |
- /* -----------------------------------------------------------------
- Programmer(s): Daniel R. Reynolds and Hilari C. Tiedeman @ SMU
- -------------------------------------------------------------------
- Copyright (c) 2011, Southern Methodist University.
- All rights reserved.
- For details, see the LICENSE file.
- -------------------------------------------------------------------
- This is the header file for the implementation of SPFGMR Krylov
- iterative linear solver. The SPFGMR algorithm is based on the
- Scaled Preconditioned Flexible GMRES (Generalized Minimal Residual)
- method [Y. Saad, SIAM J. Sci. Comput., 1993].
- The SPFGMR algorithm solves a linear system A x = b.
- Preconditioning is only allowed on the right.
- Scaling is allowed on the right, and restarts are also allowed.
- We denote the preconditioner and scaling matrices as follows:
- P = right preconditioner
- S1 = diagonal matrix of scale factors for P-inverse b
- S2 = diagonal matrix of scale factors for x
- The matrices A and P are not required explicitly; only
- routines that provide A and P-inverse as operators are required.
- In this notation, SPFGMR applies the underlying FGMRES method to
- the equivalent transformed system
- Abar xbar = bbar , where
- Abar = S1 A (P-inverse) (S2-inverse),
- bbar = S1 b , and xbar = S2 P x .
- The scaling matrix must be chosen so that the vectors S1 b and
- S2 P x have dimensionless components. If preconditioning is not
- performed (P = I), then S2 must be a scaling for x, while S1 is a
- scaling for b. Similarly, if preconditioning is performed, then S1
- must be a scaling for b, while S2 is a scaling for P x, and may
- also be taken as a scaling for b.
- The stopping test for the SPFGMR iterations is on the L2 norm of
- the scaled preconditioned residual:
- || bbar - Abar xbar ||_2 < delta
- with an input test constant delta.
- The usage of this SPFGMR solver involves supplying two routines
- and making three calls. The user-supplied routines are
- atimes (A_data, x, y) to compute y = A x, given x,
- and
- psolve (P_data, y, x, lr)
- to solve P x = y for x, given y.
- The three user calls are:
- mem = SpfgmrMalloc(lmax, vec_tmpl);
- to initialize memory,
- flag = SpfgmrSolve(mem,A_data,x,b,...,
- P_data,s1,s2,atimes,psolve,...);
- to solve the system, and
- SpfgmrFree(mem);
- to free the memory created by SpfgmrMalloc.
- Complete details for specifying atimes and psolve and for the
- usage calls are given in the paragraphs below and in iterative.h.
- -----------------------------------------------------------------*/
- #ifndef _SPFGMR_H
- #define _SPFGMR_H
- #include <sundials/sundials_iterative.h>
- #ifdef __cplusplus /* wrapper to enable C++ usage */
- extern "C" {
- #endif
- /* -----------------------------------------------------------------
- Types: SpfgmrMemRec, SpfgmrMem
- -------------------------------------------------------------------
- SpfgmrMem is a pointer to an SpfgmrMemRec which contains
- the memory needed by SpfgmrSolve. The SpfgmrMalloc routine
- returns a pointer of type SpfgmrMem which should then be passed
- in subsequent calls to SpfgmrSolve. The SpfgmrFree routine frees
- the memory allocated by SpfgmrMalloc.
- l_max is the maximum Krylov dimension that SpfgmrSolve will be
- permitted to use.
- V is the array of Krylov basis vectors v_1, ..., v_(l_max+1),
- stored in V[0], ..., V[l_max], where l_max is the second
- parameter to SpfgmrMalloc. Each v_i is a vector of type
- N_Vector.
- Z is the array of preconditioned basis vectors z_1, ...,
- z_(l_max+1), stored in Z[0], ..., Z[l_max], where l_max is the
- second parameter to SpfgmrMalloc. Each z_i is a vector of type
- N_Vector.
- Hes is the (l_max+1) x l_max Hessenberg matrix. It is stored
- row-wise so that the (i,j)th element is given by Hes[i][j].
- givens is a length 2*l_max array which represents the
- Givens rotation matrices that arise in the algorithm. The
- Givens rotation matrices F_0, F_1, ..., F_j, where F_i is
- 1
- 1
- c_i -s_i <--- row i
- s_i c_i
- 1
- 1
- are represented in the givens vector as
- givens[0]=c_0, givens[1]=s_0, givens[2]=c_1, givens[3]=s_1,
- ..., givens[2j]=c_j, givens[2j+1]=s_j.
- xcor is a vector (type N_Vector) which holds the scaled,
- preconditioned correction to the initial guess.
- yg is a length (l_max+1) array of realtype used to hold "short"
- vectors (e.g. y and g).
- vtemp is a vector (type N_Vector) used as temporary vector
- storage during calculations.
- -----------------------------------------------------------------*/
- typedef struct _SpfgmrMemRec {
- int l_max;
- N_Vector *V;
- N_Vector *Z;
- realtype **Hes;
- realtype *givens;
- N_Vector xcor;
- realtype *yg;
- N_Vector vtemp;
- } SpfgmrMemRec, *SpfgmrMem;
- /*----------------------------------------------------------------
- Function : SpfgmrMalloc
- -----------------------------------------------------------------
- SpfgmrMalloc allocates the memory used by SpfgmrSolve. It
- returns a pointer of type SpfgmrMem which the user of the
- SPGMR package should pass to SpfgmrSolve. The parameter l_max
- is the maximum Krylov dimension that SpfgmrSolve will be
- permitted to use. The parameter vec_tmpl is a pointer to an
- N_Vector used as a template to create new vectors by duplication.
- This routine returns NULL if there is a memory request failure.
- ---------------------------------------------------------------*/
- SUNDIALS_EXPORT SpfgmrMem SpfgmrMalloc(int l_max, N_Vector vec_tmpl);
- /*----------------------------------------------------------------
- Function : SpfgmrSolve
- -----------------------------------------------------------------
- SpfgmrSolve solves the linear system Ax = b using the SPFGMR
- method. The return values are given by the symbolic constants
- below. The first SpfgmrSolve parameter is a pointer to memory
- allocated by a prior call to SpfgmrMalloc.
- mem is the pointer returned by SpfgmrMalloc to the structure
- containing the memory needed by SpfgmrSolve.
- A_data is a pointer to information about the coefficient
- matrix A. This pointer is passed to the user-supplied function
- atimes.
- x is the initial guess x_0 upon entry and the solution
- N_Vector upon exit with return value SPFGMR_SUCCESS or
- SPFGMR_RES_REDUCED. For all other return values, the output x
- is undefined.
- b is the right hand side N_Vector. It is undisturbed by this
- function.
- pretype is the type of preconditioning to be used. Its
- legal possible values are enumerated in iterative.h. These
- values are PREC_NONE, PREC_LEFT, PREC_RIGHT and PREC_BOTH;
- however since this solver can only precondition on the right,
- then right-preconditioning will be done if any of the values
- PREC_LEFT, PREC_RIGHT or PREC_BOTH are provided..
- gstype is the type of Gram-Schmidt orthogonalization to be
- used. Its legal values are enumerated in iterativ.h. These
- values are MODIFIED_GS=0 and CLASSICAL_GS=1.
- delta is the tolerance on the L2 norm of the scaled,
- preconditioned residual. On return with value SPFGMR_SUCCESS,
- this residual satisfies || s1 P1_inv (b - Ax) ||_2 <= delta.
- max_restarts is the maximum number of times the algorithm is
- allowed to restart.
- maxit is the maximum number of iterations allowed within the
- solve. This value must be less than or equal to the "l_max"
- value previously supplied to SpfgmrMalloc. If maxit is too
- large, l_max will be used instead.
- P_data is a pointer to preconditioner information. This
- pointer is passed to the user-supplied function psolve.
- s1 is an N_Vector of positive scale factors for b. (Not
- tested for positivity.) Pass NULL if no scaling on b is
- required.
- s2 is an N_Vector of positive scale factors for P x, where
- P is the right preconditioner. (Not tested for positivity.)
- Pass NULL if no scaling on P x is required.
- atimes is the user-supplied function which performs the
- operation of multiplying A by a given vector. Its description
- is given in iterative.h.
- psolve is the user-supplied function which solves a
- preconditioner system Pz = r, where P is P1 or P2. Its full
- description is given in iterative.h. The psolve function will
- not be called if pretype is NONE; in that case, the user
- should pass NULL for psolve.
- res_norm is a pointer to the L2 norm of the scaled,
- preconditioned residual. On return with value SPFGMR_SUCCESS or
- SPFGMR_RES_REDUCED, (*res_norm) contains the value
- || s1 (b - Ax) ||_2 for the computed solution x.
- For all other return values, (*res_norm) is undefined. The
- caller is responsible for allocating the memory (*res_norm)
- to be filled in by SpfgmrSolve.
- nli is a pointer to the number of linear iterations done in
- the execution of SpfgmrSolve. The caller is responsible for
- allocating the memory (*nli) to be filled in by SpfgmrSolve.
- nps is a pointer to the number of calls made to psolve during
- the execution of SpfgmrSolve. The caller is responsible for
- allocating the memory (*nps) to be filled in by SpfgmrSolve.
- Note: Repeated calls can be made to SpfgmrSolve with varying
- input arguments. If, however, the problem size N or the
- maximum Krylov dimension l_max changes, then a call to
- SpfgmrMalloc must be made to obtain new memory for SpfgmrSolve
- to use.
- ---------------------------------------------------------------*/
-
- SUNDIALS_EXPORT int SpfgmrSolve(SpfgmrMem mem, void *A_data, N_Vector x,
- N_Vector b, int pretype, int gstype,
- realtype delta, int max_restarts,
- int maxit, void *P_data, N_Vector s1,
- N_Vector s2, ATimesFn atimes, PSolveFn psolve,
- realtype *res_norm, int *nli, int *nps);
- /* Return values for SpfgmrSolve */
- #define SPFGMR_SUCCESS 0 /* Converged */
- #define SPFGMR_RES_REDUCED 1 /* Did not converge, but reduced
- norm of residual */
- #define SPFGMR_CONV_FAIL 2 /* Failed to converge */
- #define SPFGMR_QRFACT_FAIL 3 /* QRfact found singular matrix */
- #define SPFGMR_PSOLVE_FAIL_REC 4 /* psolve failed recoverably */
- #define SPFGMR_ATIMES_FAIL_REC 5 /* atimes failed recoverably */
- #define SPFGMR_PSET_FAIL_REC 6 /* pset faild recoverably */
- #define SPFGMR_MEM_NULL -1 /* mem argument is NULL */
- #define SPFGMR_ATIMES_FAIL_UNREC -2 /* atimes returned failure flag */
- #define SPFGMR_PSOLVE_FAIL_UNREC -3 /* psolve failed unrecoverably */
- #define SPFGMR_GS_FAIL -4 /* Gram-Schmidt routine faiuled */
- #define SPFGMR_QRSOL_FAIL -5 /* QRsol found singular R */
- #define SPFGMR_PSET_FAIL_UNREC -6 /* pset failed unrecoverably */
- /*----------------------------------------------------------------
- Function : SpfgmrFree
- -----------------------------------------------------------------
- SpfgmrMalloc frees the memory allocated by SpfgmrMalloc. It is
- illegal to use the pointer mem after a call to SpfgmrFree.
- ---------------------------------------------------------------*/
- SUNDIALS_EXPORT void SpfgmrFree(SpfgmrMem mem);
- /*----------------------------------------------------------------
- Macro: SPFGMR_VTEMP
- -----------------------------------------------------------------
- This macro provides access to the work vector vtemp in the
- memory block of the SPFGMR module. The argument mem is the
- memory pointer returned by SpfgmrMalloc, of type SpfgmrMem,
- and the macro value is of type N_Vector.
- On a return from SpfgmrSolve with *nli = 0, this vector
- contains the scaled preconditioned initial residual,
- s1 * P1_inverse * (b - A x_0).
- ---------------------------------------------------------------*/
- #define SPFGMR_VTEMP(mem) (mem->vtemp)
- #ifdef __cplusplus
- }
- #endif
- #endif
|