[16222] | 1 | /* |
---|
| 2 | * ----------------------------------------------------------------- |
---|
| 3 | * $Revision$ |
---|
| 4 | * $Date$ |
---|
| 5 | * ----------------------------------------------------------------- |
---|
| 6 | * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and |
---|
| 7 | * Radu Serban @ LLNL |
---|
| 8 | * ----------------------------------------------------------------- |
---|
| 9 | * LLNS Copyright Start |
---|
| 10 | * Copyright (c) 2014, Lawrence Livermore National Security |
---|
| 11 | * This work was performed under the auspices of the U.S. Department |
---|
| 12 | * of Energy by Lawrence Livermore National Laboratory in part under |
---|
| 13 | * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344. |
---|
| 14 | * Produced at the Lawrence Livermore National Laboratory. |
---|
| 15 | * All rights reserved. |
---|
| 16 | * For details, see the LICENSE file. |
---|
| 17 | * LLNS Copyright End |
---|
| 18 | * ----------------------------------------------------------------- |
---|
| 19 | * This is the header file for the implementation of SPGMR Krylov |
---|
| 20 | * iterative linear solver. The SPGMR algorithm is based on the |
---|
| 21 | * Scaled Preconditioned GMRES (Generalized Minimal Residual) |
---|
| 22 | * method. |
---|
| 23 | * |
---|
| 24 | * The SPGMR algorithm solves a linear system A x = b. |
---|
| 25 | * Preconditioning is allowed on the left, right, or both. |
---|
| 26 | * Scaling is allowed on both sides, and restarts are also allowed. |
---|
| 27 | * We denote the preconditioner and scaling matrices as follows: |
---|
| 28 | * P1 = left preconditioner |
---|
| 29 | * P2 = right preconditioner |
---|
| 30 | * S1 = diagonal matrix of scale factors for P1-inverse b |
---|
| 31 | * S2 = diagonal matrix of scale factors for P2 x |
---|
| 32 | * The matrices A, P1, and P2 are not required explicitly; only |
---|
| 33 | * routines that provide A, P1-inverse, and P2-inverse as |
---|
| 34 | * operators are required. |
---|
| 35 | * |
---|
| 36 | * In this notation, SPGMR applies the underlying GMRES method to |
---|
| 37 | * the equivalent transformed system |
---|
| 38 | * Abar xbar = bbar , where |
---|
| 39 | * Abar = S1 (P1-inverse) A (P2-inverse) (S2-inverse) , |
---|
| 40 | * bbar = S1 (P1-inverse) b , and xbar = S2 P2 x . |
---|
| 41 | * |
---|
| 42 | * The scaling matrices must be chosen so that vectors S1 |
---|
| 43 | * P1-inverse b and S2 P2 x have dimensionless components. |
---|
| 44 | * If preconditioning is done on the left only (P2 = I), by a |
---|
| 45 | * matrix P, then S2 must be a scaling for x, while S1 is a |
---|
| 46 | * scaling for P-inverse b, and so may also be taken as a scaling |
---|
| 47 | * for x. Similarly, if preconditioning is done on the right only |
---|
| 48 | * (P1 = I, P2 = P), then S1 must be a scaling for b, while S2 is |
---|
| 49 | * a scaling for P x, and may also be taken as a scaling for b. |
---|
| 50 | * |
---|
| 51 | * The stopping test for the SPGMR iterations is on the L2 norm of |
---|
| 52 | * the scaled preconditioned residual: |
---|
| 53 | * || bbar - Abar xbar ||_2 < delta |
---|
| 54 | * with an input test constant delta. |
---|
| 55 | * |
---|
| 56 | * The usage of this SPGMR solver involves supplying two routines |
---|
| 57 | * and making three calls. The user-supplied routines are |
---|
| 58 | * atimes (A_data, x, y) to compute y = A x, given x, |
---|
| 59 | * and |
---|
| 60 | * psolve (P_data, y, x, lr) |
---|
| 61 | * to solve P1 x = y or P2 x = y for x, given y. |
---|
| 62 | * The three user calls are: |
---|
| 63 | * mem = SpgmrMalloc(lmax, vec_tmpl); |
---|
| 64 | * to initialize memory, |
---|
| 65 | * flag = SpgmrSolve(mem,A_data,x,b,..., |
---|
| 66 | * P_data,s1,s2,atimes,psolve,...); |
---|
| 67 | * to solve the system, and |
---|
| 68 | * SpgmrFree(mem); |
---|
| 69 | * to free the memory created by SpgmrMalloc. |
---|
| 70 | * Complete details for specifying atimes and psolve and for the |
---|
| 71 | * usage calls are given below and in sundials_iterative.h. |
---|
| 72 | * ----------------------------------------------------------------- |
---|
| 73 | */ |
---|
| 74 | |
---|
| 75 | #ifndef _SPGMR_H |
---|
| 76 | #define _SPGMR_H |
---|
| 77 | |
---|
| 78 | #include <sundials/sundials_iterative.h> |
---|
| 79 | |
---|
| 80 | #ifdef __cplusplus /* wrapper to enable C++ usage */ |
---|
| 81 | extern "C" { |
---|
| 82 | #endif |
---|
| 83 | |
---|
| 84 | /* |
---|
| 85 | * ----------------------------------------------------------------- |
---|
| 86 | * Types: SpgmrMemRec, SpgmrMem |
---|
| 87 | * ----------------------------------------------------------------- |
---|
| 88 | * SpgmrMem is a pointer to an SpgmrMemRec which contains |
---|
| 89 | * the memory needed by SpgmrSolve. The SpgmrMalloc routine |
---|
| 90 | * returns a pointer of type SpgmrMem which should then be passed |
---|
| 91 | * in subsequent calls to SpgmrSolve. The SpgmrFree routine frees |
---|
| 92 | * the memory allocated by SpgmrMalloc. |
---|
| 93 | * |
---|
| 94 | * l_max is the maximum Krylov dimension that SpgmrSolve will be |
---|
| 95 | * permitted to use. |
---|
| 96 | * |
---|
| 97 | * V is the array of Krylov basis vectors v_1, ..., v_(l_max+1), |
---|
| 98 | * stored in V[0], ..., V[l_max], where l_max is the second |
---|
| 99 | * parameter to SpgmrMalloc. Each v_i is a vector of type |
---|
| 100 | * N_Vector. |
---|
| 101 | * |
---|
| 102 | * Hes is the (l_max+1) x l_max Hessenberg matrix. It is stored |
---|
| 103 | * row-wise so that the (i,j)th element is given by Hes[i][j]. |
---|
| 104 | * |
---|
| 105 | * givens is a length 2*l_max array which represents the |
---|
| 106 | * Givens rotation matrices that arise in the algorithm. The |
---|
| 107 | * Givens rotation matrices F_0, F_1, ..., F_j, where F_i is |
---|
| 108 | * |
---|
| 109 | * 1 |
---|
| 110 | * 1 |
---|
| 111 | * c_i -s_i <--- row i |
---|
| 112 | * s_i c_i |
---|
| 113 | * 1 |
---|
| 114 | * 1 |
---|
| 115 | * |
---|
| 116 | * are represented in the givens vector as |
---|
| 117 | * givens[0]=c_0, givens[1]=s_0, givens[2]=c_1, givens[3]=s_1, |
---|
| 118 | * ..., givens[2j]=c_j, givens[2j+1]=s_j. |
---|
| 119 | * |
---|
| 120 | * xcor is a vector (type N_Vector) which holds the scaled, |
---|
| 121 | * preconditioned correction to the initial guess. |
---|
| 122 | * |
---|
| 123 | * yg is a length (l_max+1) array of realtype used to hold "short" |
---|
| 124 | * vectors (e.g. y and g). |
---|
| 125 | * |
---|
| 126 | * vtemp is a vector (type N_Vector) used as temporary vector |
---|
| 127 | * storage during calculations. |
---|
| 128 | * ----------------------------------------------------------------- |
---|
| 129 | */ |
---|
| 130 | |
---|
| 131 | typedef struct _SpgmrMemRec { |
---|
| 132 | |
---|
| 133 | int l_max; |
---|
| 134 | |
---|
| 135 | N_Vector *V; |
---|
| 136 | realtype **Hes; |
---|
| 137 | realtype *givens; |
---|
| 138 | N_Vector xcor; |
---|
| 139 | realtype *yg; |
---|
| 140 | N_Vector vtemp; |
---|
| 141 | |
---|
| 142 | } SpgmrMemRec, *SpgmrMem; |
---|
| 143 | |
---|
| 144 | /* |
---|
| 145 | * ----------------------------------------------------------------- |
---|
| 146 | * Function : SpgmrMalloc |
---|
| 147 | * ----------------------------------------------------------------- |
---|
| 148 | * SpgmrMalloc allocates the memory used by SpgmrSolve. It |
---|
| 149 | * returns a pointer of type SpgmrMem which the user of the |
---|
| 150 | * SPGMR package should pass to SpgmrSolve. The parameter l_max |
---|
| 151 | * is the maximum Krylov dimension that SpgmrSolve will be |
---|
| 152 | * permitted to use. The parameter vec_tmpl is a pointer to an |
---|
| 153 | * N_Vector used as a template to create new vectors by duplication. |
---|
| 154 | * This routine returns NULL if there is a memory request failure. |
---|
| 155 | * ----------------------------------------------------------------- |
---|
| 156 | */ |
---|
| 157 | |
---|
| 158 | SUNDIALS_EXPORT SpgmrMem SpgmrMalloc(int l_max, N_Vector vec_tmpl); |
---|
| 159 | |
---|
| 160 | /* |
---|
| 161 | * ----------------------------------------------------------------- |
---|
| 162 | * Function : SpgmrSolve |
---|
| 163 | * ----------------------------------------------------------------- |
---|
| 164 | * SpgmrSolve solves the linear system Ax = b using the SPGMR |
---|
| 165 | * method. The return values are given by the symbolic constants |
---|
| 166 | * below. The first SpgmrSolve parameter is a pointer to memory |
---|
| 167 | * allocated by a prior call to SpgmrMalloc. |
---|
| 168 | * |
---|
| 169 | * mem is the pointer returned by SpgmrMalloc to the structure |
---|
| 170 | * containing the memory needed by SpgmrSolve. |
---|
| 171 | * |
---|
| 172 | * A_data is a pointer to information about the coefficient |
---|
| 173 | * matrix A. This pointer is passed to the user-supplied function |
---|
| 174 | * atimes. |
---|
| 175 | * |
---|
| 176 | * x is the initial guess x_0 upon entry and the solution |
---|
| 177 | * N_Vector upon exit with return value SPGMR_SUCCESS or |
---|
| 178 | * SPGMR_RES_REDUCED. For all other return values, the output x |
---|
| 179 | * is undefined. |
---|
| 180 | * |
---|
| 181 | * b is the right hand side N_Vector. It is undisturbed by this |
---|
| 182 | * function. |
---|
| 183 | * |
---|
| 184 | * pretype is the type of preconditioning to be used. Its |
---|
| 185 | * legal values are enumerated in sundials_iterative.h. These |
---|
| 186 | * values are PREC_NONE=0, PREC_LEFT=1, PREC_RIGHT=2, and |
---|
| 187 | * PREC_BOTH=3. |
---|
| 188 | * |
---|
| 189 | * gstype is the type of Gram-Schmidt orthogonalization to be |
---|
| 190 | * used. Its legal values are enumerated in sundials_iterative.h. |
---|
| 191 | * These values are MODIFIED_GS=0 and CLASSICAL_GS=1. |
---|
| 192 | * |
---|
| 193 | * delta is the tolerance on the L2 norm of the scaled, |
---|
| 194 | * preconditioned residual. On return with value SPGMR_SUCCESS, |
---|
| 195 | * this residual satisfies || s1 P1_inv (b - Ax) ||_2 <= delta. |
---|
| 196 | * |
---|
| 197 | * max_restarts is the maximum number of times the algorithm is |
---|
| 198 | * allowed to restart. |
---|
| 199 | * |
---|
| 200 | * P_data is a pointer to preconditioner information. This |
---|
| 201 | * pointer is passed to the user-supplied function psolve. |
---|
| 202 | * |
---|
| 203 | * s1 is an N_Vector of positive scale factors for P1-inv b, where |
---|
| 204 | * P1 is the left preconditioner. (Not tested for positivity.) |
---|
| 205 | * Pass NULL if no scaling on P1-inv b is required. |
---|
| 206 | * |
---|
| 207 | * s2 is an N_Vector of positive scale factors for P2 x, where |
---|
| 208 | * P2 is the right preconditioner. (Not tested for positivity.) |
---|
| 209 | * Pass NULL if no scaling on P2 x is required. |
---|
| 210 | * |
---|
| 211 | * atimes is the user-supplied function which performs the |
---|
| 212 | * operation of multiplying A by a given vector. Its description |
---|
| 213 | * is given in sundials_iterative.h. |
---|
| 214 | * |
---|
| 215 | * psolve is the user-supplied function which solves a |
---|
| 216 | * preconditioner system Pz = r, where P is P1 or P2. Its full |
---|
| 217 | * description is given in sundials_iterative.h. The psolve function |
---|
| 218 | * will not be called if pretype is NONE; in that case, the user |
---|
| 219 | * should pass NULL for psolve. |
---|
| 220 | * |
---|
| 221 | * res_norm is a pointer to the L2 norm of the scaled, |
---|
| 222 | * preconditioned residual. On return with value SPGMR_SUCCESS or |
---|
| 223 | * SPGMR_RES_REDUCED, (*res_norm) contains the value |
---|
| 224 | * || s1 P1_inv (b - Ax) ||_2 for the computed solution x. |
---|
| 225 | * For all other return values, (*res_norm) is undefined. The |
---|
| 226 | * caller is responsible for allocating the memory (*res_norm) |
---|
| 227 | * to be filled in by SpgmrSolve. |
---|
| 228 | * |
---|
| 229 | * nli is a pointer to the number of linear iterations done in |
---|
| 230 | * the execution of SpgmrSolve. The caller is responsible for |
---|
| 231 | * allocating the memory (*nli) to be filled in by SpgmrSolve. |
---|
| 232 | * |
---|
| 233 | * nps is a pointer to the number of calls made to psolve during |
---|
| 234 | * the execution of SpgmrSolve. The caller is responsible for |
---|
| 235 | * allocating the memory (*nps) to be filled in by SpgmrSolve. |
---|
| 236 | * |
---|
| 237 | * Note: Repeated calls can be made to SpgmrSolve with varying |
---|
| 238 | * input arguments. If, however, the problem size N or the |
---|
| 239 | * maximum Krylov dimension l_max changes, then a call to |
---|
| 240 | * SpgmrMalloc must be made to obtain new memory for SpgmrSolve |
---|
| 241 | * to use. |
---|
| 242 | * ----------------------------------------------------------------- |
---|
| 243 | */ |
---|
| 244 | |
---|
| 245 | SUNDIALS_EXPORT int SpgmrSolve(SpgmrMem mem, void *A_data, N_Vector x, N_Vector b, |
---|
| 246 | int pretype, int gstype, realtype delta, |
---|
| 247 | int max_restarts, void *P_data, N_Vector s1, |
---|
| 248 | N_Vector s2, ATimesFn atimes, PSolveFn psolve, |
---|
| 249 | realtype *res_norm, int *nli, int *nps); |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | /* Return values for SpgmrSolve */ |
---|
| 253 | |
---|
| 254 | #define SPGMR_SUCCESS 0 /* Converged */ |
---|
| 255 | #define SPGMR_RES_REDUCED 1 /* Did not converge, but reduced |
---|
| 256 | norm of residual */ |
---|
| 257 | #define SPGMR_CONV_FAIL 2 /* Failed to converge */ |
---|
| 258 | #define SPGMR_QRFACT_FAIL 3 /* QRfact found singular matrix */ |
---|
| 259 | #define SPGMR_PSOLVE_FAIL_REC 4 /* psolve failed recoverably */ |
---|
| 260 | #define SPGMR_ATIMES_FAIL_REC 5 /* atimes failed recoverably */ |
---|
| 261 | #define SPGMR_PSET_FAIL_REC 6 /* pset faild recoverably */ |
---|
| 262 | |
---|
| 263 | #define SPGMR_MEM_NULL -1 /* mem argument is NULL */ |
---|
| 264 | #define SPGMR_ATIMES_FAIL_UNREC -2 /* atimes returned failure flag */ |
---|
| 265 | #define SPGMR_PSOLVE_FAIL_UNREC -3 /* psolve failed unrecoverably */ |
---|
| 266 | #define SPGMR_GS_FAIL -4 /* Gram-Schmidt routine faiuled */ |
---|
| 267 | #define SPGMR_QRSOL_FAIL -5 /* QRsol found singular R */ |
---|
| 268 | #define SPGMR_PSET_FAIL_UNREC -6 /* pset failed unrecoverably */ |
---|
| 269 | |
---|
| 270 | /* |
---|
| 271 | * ----------------------------------------------------------------- |
---|
| 272 | * Function : SpgmrFree |
---|
| 273 | * ----------------------------------------------------------------- |
---|
| 274 | * SpgmrMalloc frees the memory allocated by SpgmrMalloc. It is |
---|
| 275 | * illegal to use the pointer mem after a call to SpgmrFree. |
---|
| 276 | * ----------------------------------------------------------------- |
---|
| 277 | */ |
---|
| 278 | |
---|
| 279 | SUNDIALS_EXPORT void SpgmrFree(SpgmrMem mem); |
---|
| 280 | |
---|
| 281 | /* |
---|
| 282 | * ----------------------------------------------------------------- |
---|
| 283 | * Macro: SPGMR_VTEMP |
---|
| 284 | * ----------------------------------------------------------------- |
---|
| 285 | * This macro provides access to the work vector vtemp in the |
---|
| 286 | * memory block of the SPGMR module. The argument mem is the |
---|
| 287 | * memory pointer returned by SpgmrMalloc, of type SpgmrMem, |
---|
| 288 | * and the macro value is of type N_Vector. |
---|
| 289 | * On a return from SpgmrSolve with *nli = 0, this vector |
---|
| 290 | * contains the scaled preconditioned initial residual, |
---|
| 291 | * s1 * P1_inverse * (b - A x_0). |
---|
| 292 | * ----------------------------------------------------------------- |
---|
| 293 | */ |
---|
| 294 | |
---|
| 295 | #define SPGMR_VTEMP(mem) (mem->vtemp) |
---|
| 296 | |
---|
| 297 | #ifdef __cplusplus |
---|
| 298 | } |
---|
| 299 | #endif |
---|
| 300 | |
---|
| 301 | #endif |
---|