1 /* 2 * ----------------------------------------------------------------- 3 * $Revision: 4378 $ 4 * $Date: 2015-02-19 10:55:14 -0800 (Thu, 19 Feb 2015) $ 5 * ----------------------------------------------------------------- 6 * Programmer(s): Radu Serban @ LLNL 7 * ----------------------------------------------------------------- 8 * LLNS Copyright Start 9 * Copyright (c) 2014, Lawrence Livermore National Security 10 * This work was performed under the auspices of the U.S. Department 11 * of Energy by Lawrence Livermore National Laboratory in part under 12 * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344. 13 * Produced at the Lawrence Livermore National Laboratory. 14 * All rights reserved. 15 * For details, see the LICENSE file. 16 * LLNS Copyright End 17 * ----------------------------------------------------------------- 18 * This is the common header file for the Scaled, Preconditioned 19 * Iterative Linear Solvers in CVODES. 20 * 21 * Part I contains type definitions and functions for using the 22 * iterative linear solvers on forward problems 23 * (IVP integration and/or FSA) 24 * 25 * Part II contains type definitions and functions for using the 26 * iterative linear solvers on adjoint (backward) problems 27 * ----------------------------------------------------------------- 28 */ 29 30 #ifndef _CVSSPILS_H 31 #define _CVSSPILS_H 32 33 #include <sundials/sundials_iterative.h> 34 #include <sundials/sundials_nvector.h> 35 36 #ifdef __cplusplus /* wrapper to enable C++ usage */ 37 extern "C" { 38 #endif 39 40 /* 41 * ----------------------------------------------------------------- 42 * CVSPILS return values 43 * ----------------------------------------------------------------- 44 */ 45 46 #define CVSPILS_SUCCESS 0 47 #define CVSPILS_MEM_NULL -1 48 #define CVSPILS_LMEM_NULL -2 49 #define CVSPILS_ILL_INPUT -3 50 #define CVSPILS_MEM_FAIL -4 51 #define CVSPILS_PMEM_NULL -5 52 53 /* Return values for the adjoint module */ 54 55 #define CVSPILS_NO_ADJ -101 56 #define CVSPILS_LMEMB_NULL -102 57 58 /* 59 * ----------------------------------------------------------------- 60 * CVSPILS solver constants 61 * ----------------------------------------------------------------- 62 * CVSPILS_MAXL : default value for the maximum Krylov 63 * dimension 64 * 65 * CVSPILS_MSBPRE : maximum number of steps between 66 * preconditioner evaluations 67 * 68 * CVSPILS_DGMAX : maximum change in gamma between 69 * preconditioner evaluations 70 * 71 * CVSPILS_EPLIN : default value for factor by which the 72 * tolerance on the nonlinear iteration is 73 * multiplied to get a tolerance on the linear 74 * iteration 75 * ----------------------------------------------------------------- 76 */ 77 78 #define CVSPILS_MAXL 5 79 #define CVSPILS_MSBPRE 50 80 #define CVSPILS_DGMAX RCONST(0.2) 81 #define CVSPILS_EPLIN RCONST(0.05) 82 83 /* 84 * ----------------------------------------------------------------- 85 * PART I - forward problems 86 * ----------------------------------------------------------------- 87 */ 88 89 /* 90 * ----------------------------------------------------------------- 91 * Type : CVSpilsPrecSetupFn 92 * ----------------------------------------------------------------- 93 * The user-supplied preconditioner setup function PrecSetup and 94 * the user-supplied preconditioner solve function PrecSolve 95 * together must define left and right preconditoner matrices 96 * P1 and P2 (either of which may be trivial), such that the 97 * product P1*P2 is an approximation to the Newton matrix 98 * M = I - gamma*J. Here J is the system Jacobian J = df/dy, 99 * and gamma is a scalar proportional to the integration step 100 * size h. The solution of systems P z = r, with P = P1 or P2, 101 * is to be carried out by the PrecSolve function, and PrecSetup 102 * is to do any necessary setup operations. 103 * 104 * The user-supplied preconditioner setup function PrecSetup 105 * is to evaluate and preprocess any Jacobian-related data 106 * needed by the preconditioner solve function PrecSolve. 107 * This might include forming a crude approximate Jacobian, 108 * and performing an LU factorization on the resulting 109 * approximation to M. This function will not be called in 110 * advance of every call to PrecSolve, but instead will be called 111 * only as often as necessary to achieve convergence within the 112 * Newton iteration. If the PrecSolve function needs no 113 * preparation, the PrecSetup function can be NULL. 114 * 115 * For greater efficiency, the PrecSetup function may save 116 * Jacobian-related data and reuse it, rather than generating it 117 * from scratch. In this case, it should use the input flag jok 118 * to decide whether to recompute the data, and set the output 119 * flag *jcurPtr accordingly. 120 * 121 * Each call to the PrecSetup function is preceded by a call to 122 * the RhsFn f with the same (t,y) arguments. Thus the PrecSetup 123 * function can use any auxiliary data that is computed and 124 * saved by the f function and made accessible to PrecSetup. 125 * 126 * A function PrecSetup must have the prototype given below. 127 * Its parameters are as follows: 128 * 129 * t is the current value of the independent variable. 130 * 131 * y is the current value of the dependent variable vector, 132 * namely the predicted value of y(t). 133 * 134 * fy is the vector f(t,y). 135 * 136 * jok is an input flag indicating whether Jacobian-related 137 * data needs to be recomputed, as follows: 138 * jok == FALSE means recompute Jacobian-related data 139 * from scratch. 140 * jok == TRUE means that Jacobian data, if saved from 141 * the previous PrecSetup call, can be reused 142 * (with the current value of gamma). 143 * A Precset call with jok == TRUE can only occur after 144 * a call with jok == FALSE. 145 * 146 * jcurPtr is a pointer to an output integer flag which is 147 * to be set by PrecSetup as follows: 148 * Set *jcurPtr = TRUE if Jacobian data was recomputed. 149 * Set *jcurPtr = FALSE if Jacobian data was not recomputed, 150 * but saved data was reused. 151 * 152 * gamma is the scalar appearing in the Newton matrix. 153 * 154 * user_data is a pointer to user data - the same as the user_data 155 * parameter passed to the CVodeSetUserData function. 156 * 157 * tmp1, tmp2, and tmp3 are pointers to memory allocated 158 * for N_Vectors which can be used by 159 * CVSpilsPrecSetupFn as temporary storage or 160 * work space. 161 * 162 * NOTE: If the user's preconditioner needs other quantities, 163 * they are accessible as follows: hcur (the current stepsize) 164 * and ewt (the error weight vector) are accessible through 165 * CVodeGetCurrentStep and CVodeGetErrWeights, respectively). 166 * The unit roundoff is available as UNIT_ROUNDOFF defined in 167 * sundials_types.h. 168 * 169 * Returned value: 170 * The value to be returned by the PrecSetup function is a flag 171 * indicating whether it was successful. This value should be 172 * 0 if successful, 173 * > 0 for a recoverable error (step will be retried), 174 * < 0 for an unrecoverable error (integration is halted). 175 * ----------------------------------------------------------------- 176 */ 177 178 typedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy, 179 booleantype jok, booleantype *jcurPtr, 180 realtype gamma, void *user_data, 181 N_Vector tmp1, N_Vector tmp2, 182 N_Vector tmp3); 183 184 /* 185 * ----------------------------------------------------------------- 186 * Type : CVSpilsPrecSolveFn 187 * ----------------------------------------------------------------- 188 * The user-supplied preconditioner solve function PrecSolve 189 * is to solve a linear system P z = r in which the matrix P is 190 * one of the preconditioner matrices P1 or P2, depending on the 191 * type of preconditioning chosen. 192 * 193 * A function PrecSolve must have the prototype given below. 194 * Its parameters are as follows: 195 * 196 * t is the current value of the independent variable. 197 * 198 * y is the current value of the dependent variable vector. 199 * 200 * fy is the vector f(t,y). 201 * 202 * r is the right-hand side vector of the linear system. 203 * 204 * z is the output vector computed by PrecSolve. 205 * 206 * gamma is the scalar appearing in the Newton matrix. 207 * 208 * delta is an input tolerance for use by PSolve if it uses 209 * an iterative method in its solution. In that case, 210 * the residual vector Res = r - P z of the system 211 * should be made less than delta in weighted L2 norm, 212 * i.e., sqrt [ Sum (Res[i]*ewt[i])^2 ] < delta. 213 * Note: the error weight vector ewt can be obtained 214 * through a call to the routine CVodeGetErrWeights. 215 * 216 * lr is an input flag indicating whether PrecSolve is to use 217 * the left preconditioner P1 or right preconditioner 218 * P2: lr = 1 means use P1, and lr = 2 means use P2. 219 * 220 * user_data is a pointer to user data - the same as the user_data 221 * parameter passed to the CVodeSetUserData function. 222 * 223 * tmp is a pointer to memory allocated for an N_Vector 224 * which can be used by PSolve for work space. 225 * 226 * Returned value: 227 * The value to be returned by the PrecSolve function is a flag 228 * indicating whether it was successful. This value should be 229 * 0 if successful, 230 * positive for a recoverable error (step will be retried), 231 * negative for an unrecoverable error (integration is halted). 232 * ----------------------------------------------------------------- 233 */ 234 235 typedef int (*CVSpilsPrecSolveFn)(realtype t, N_Vector y, N_Vector fy, 236 N_Vector r, N_Vector z, 237 realtype gamma, realtype delta, 238 int lr, void *user_data, N_Vector tmp); 239 240 /* 241 * ----------------------------------------------------------------- 242 * Type : CVSpilsJacTimesVecFn 243 * ----------------------------------------------------------------- 244 * The user-supplied function jtimes is to generate the product 245 * J*v for given v, where J is the Jacobian df/dy, or an 246 * approximation to it, and v is a given vector. It should return 247 * 0 if successful a positive value for a recoverable error or 248 * a negative value for an unrecoverable failure. 249 * 250 * A function jtimes must have the prototype given below. Its 251 * parameters are as follows: 252 * 253 * v is the N_Vector to be multiplied by J. 254 * 255 * Jv is the output N_Vector containing J*v. 256 * 257 * t is the current value of the independent variable. 258 * 259 * y is the current value of the dependent variable 260 * vector. 261 * 262 * fy is the vector f(t,y). 263 * 264 * user_data is a pointer to user data, the same as the user_data 265 * parameter passed to the CVodeSetUserData function. 266 * 267 * tmp is a pointer to memory allocated for an N_Vector 268 * which can be used by Jtimes for work space. 269 * ----------------------------------------------------------------- 270 */ 271 272 typedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t, 273 N_Vector y, N_Vector fy, 274 void *user_data, N_Vector tmp); 275 276 277 /* 278 * ----------------------------------------------------------------- 279 * Optional inputs to the CVSPILS linear solver 280 * ----------------------------------------------------------------- 281 * 282 * CVSpilsSetPrecType resets the type of preconditioner, pretype, 283 * from the value previously set. 284 * This must be one of PREC_NONE, PREC_LEFT, 285 * PREC_RIGHT, or PREC_BOTH. 286 * 287 * CVSpilsSetGSType specifies the type of Gram-Schmidt 288 * orthogonalization to be used. This must be one of 289 * the two enumeration constants MODIFIED_GS or 290 * CLASSICAL_GS defined in iterative.h. These correspond 291 * to using modified Gram-Schmidt and classical 292 * Gram-Schmidt, respectively. 293 * Default value is MODIFIED_GS. 294 * 295 * CVSpilsSetMaxl resets the maximum Krylov subspace size, maxl, 296 * from the value previously set. 297 * An input value <= 0, gives the default value. 298 * 299 * CVSpilsSetEpsLin specifies the factor by which the tolerance on 300 * the nonlinear iteration is multiplied to get a 301 * tolerance on the linear iteration. 302 * Default value is 0.05. 303 * 304 * CVSpilsSetPreconditioner specifies the PrecSetup and PrecSolve functions. 305 * Default is NULL for both arguments (no preconditioning). 306 * 307 * CVSpilsSetJacTimesVecFn specifies the jtimes function. Default is to use 308 * an internal finite difference approximation routine. 309 * 310 * The return value of CVSpilsSet* is one of: 311 * CVSPILS_SUCCESS if successful 312 * CVSPILS_MEM_NULL if the cvode memory was NULL 313 * CVSPILS_LMEM_NULL if the linear solver memory was NULL 314 * CVSPILS_ILL_INPUT if an input has an illegal value 315 * ----------------------------------------------------------------- 316 */ 317 318 SUNDIALS_EXPORT int CVSpilsSetPrecType(void *cvode_mem, int pretype); 319 SUNDIALS_EXPORT int CVSpilsSetGSType(void *cvode_mem, int gstype); 320 SUNDIALS_EXPORT int CVSpilsSetMaxl(void *cvode_mem, int maxl); 321 SUNDIALS_EXPORT int CVSpilsSetEpsLin(void *cvode_mem, realtype eplifac); 322 SUNDIALS_EXPORT int CVSpilsSetPreconditioner(void *cvode_mem, 323 CVSpilsPrecSetupFn pset, 324 CVSpilsPrecSolveFn psolve); 325 SUNDIALS_EXPORT int CVSpilsSetJacTimesVecFn(void *cvode_mem, 326 CVSpilsJacTimesVecFn jtv); 327 328 /* 329 * ----------------------------------------------------------------- 330 * Optional outputs from the CVSPILS linear solver 331 * ----------------------------------------------------------------- 332 * CVSpilsGetWorkSpace returns the real and integer workspace used 333 * by the SPILS module. 334 * 335 * CVSpilsGetNumPrecEvals returns the number of preconditioner 336 * evaluations, i.e. the number of calls made 337 * to PrecSetup with jok==FALSE. 338 * 339 * CVSpilsGetNumPrecSolves returns the number of calls made to 340 * PrecSolve. 341 * 342 * CVSpilsGetNumLinIters returns the number of linear iterations. 343 * 344 * CVSpilsGetNumConvFails returns the number of linear 345 * convergence failures. 346 * 347 * CVSpilsGetNumJtimesEvals returns the number of calls to jtimes. 348 * 349 * CVSpilsGetNumRhsEvals returns the number of calls to the user 350 * f routine due to finite difference Jacobian 351 * times vector evaluation. 352 * 353 * CVSpilsGetLastFlag returns the last error flag set by any of 354 * the CVSPILS interface functions. 355 * 356 * The return value of CVSpilsGet* is one of: 357 * CVSPILS_SUCCESS if successful 358 * CVSPILS_MEM_NULL if the cvode memory was NULL 359 * CVSPILS_LMEM_NULL if the linear solver memory was NULL 360 * ----------------------------------------------------------------- 361 */ 362 363 SUNDIALS_EXPORT int CVSpilsGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS); 364 SUNDIALS_EXPORT int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals); 365 SUNDIALS_EXPORT int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves); 366 SUNDIALS_EXPORT int CVSpilsGetNumLinIters(void *cvode_mem, long int *nliters); 367 SUNDIALS_EXPORT int CVSpilsGetNumConvFails(void *cvode_mem, long int *nlcfails); 368 SUNDIALS_EXPORT int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals); 369 SUNDIALS_EXPORT int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS); 370 SUNDIALS_EXPORT int CVSpilsGetLastFlag(void *cvode_mem, long int *flag); 371 372 /* 373 * ----------------------------------------------------------------- 374 * The following function returns the name of the constant 375 * associated with a CVSPILS return flag 376 * ----------------------------------------------------------------- 377 */ 378 379 SUNDIALS_EXPORT char *CVSpilsGetReturnFlagName(long int flag); 380 381 382 /* 383 * ----------------------------------------------------------------- 384 * PART II - backward problems 385 * ----------------------------------------------------------------- 386 */ 387 388 /* 389 * ----------------------------------------------------------------- 390 * Type : CVSpilsPrecSetupFnB 391 * ----------------------------------------------------------------- 392 * A function PrecSetupB for the adjoint (backward) problem must have 393 * the prototype given below. 394 * ----------------------------------------------------------------- 395 */ 396 397 typedef int (*CVSpilsPrecSetupFnB)(realtype t, N_Vector y, 398 N_Vector yB, N_Vector fyB, 399 booleantype jokB, 400 booleantype *jcurPtrB, realtype gammaB, 401 void *user_dataB, 402 N_Vector tmp1B, N_Vector tmp2B, 403 N_Vector tmp3B); 404 405 406 /* 407 * ----------------------------------------------------------------- 408 * Type : CVSpilsPrecSetupFnBS 409 * ----------------------------------------------------------------- 410 * A function PrecSetupBS for the adjoint (backward) problem must have 411 * the prototype given below. 412 * ----------------------------------------------------------------- 413 */ 414 415 typedef int (*CVSpilsPrecSetupFnBS)(realtype t, N_Vector y, N_Vector *yS, 416 N_Vector yB, N_Vector fyB, 417 booleantype jokB, 418 booleantype *jcurPtrB, realtype gammaB, 419 void *user_dataB, 420 N_Vector tmp1B, N_Vector tmp2B, 421 N_Vector tmp3B); 422 423 424 /* 425 * ----------------------------------------------------------------- 426 * Type : CVSpilsPrecSolveFnB 427 * ----------------------------------------------------------------- 428 * A function PrecSolveB for the adjoint (backward) problem must 429 * have the prototype given below. 430 * ----------------------------------------------------------------- 431 */ 432 433 typedef int (*CVSpilsPrecSolveFnB)(realtype t, N_Vector y, 434 N_Vector yB, N_Vector fyB, 435 N_Vector rB, N_Vector zB, 436 realtype gammaB, realtype deltaB, 437 int lrB, void *user_dataB, N_Vector tmpB); 438 439 /* 440 * ----------------------------------------------------------------- 441 * Type : CVSpilsPrecSolveFnBS 442 * ----------------------------------------------------------------- 443 * A function PrecSolveBS for the adjoint (backward) problem must 444 * have the prototype given below. 445 * ----------------------------------------------------------------- 446 */ 447 448 typedef int (*CVSpilsPrecSolveFnBS)(realtype t, N_Vector y, N_Vector *yS, 449 N_Vector yB, N_Vector fyB, 450 N_Vector rB, N_Vector zB, 451 realtype gammaB, realtype deltaB, 452 int lrB, void *user_dataB, N_Vector tmpB); 453 454 /* 455 * ----------------------------------------------------------------- 456 * Type : CVSpilsJacTimesVecFnB 457 * ----------------------------------------------------------------- 458 * A function jtimesB for the adjoint (backward) problem must have 459 * the prototype given below. 460 * ----------------------------------------------------------------- 461 */ 462 463 typedef int (*CVSpilsJacTimesVecFnB)(N_Vector vB, N_Vector JvB, realtype t, 464 N_Vector y, N_Vector yB, N_Vector fyB, 465 void *jac_dataB, N_Vector tmpB); 466 467 /* 468 * ----------------------------------------------------------------- 469 * Type : CVSpilsJacTimesVecFnBS 470 * ----------------------------------------------------------------- 471 * A function jtimesBS for the adjoint (backward) problem must have 472 * the prototype given below. 473 * ----------------------------------------------------------------- 474 */ 475 476 typedef int (*CVSpilsJacTimesVecFnBS)(N_Vector vB, N_Vector JvB, 477 realtype t, N_Vector y, N_Vector *yS, 478 N_Vector yB, N_Vector fyB, 479 void *jac_dataB, N_Vector tmpB); 480 481 /* 482 * ----------------------------------------------------------------- 483 * Functions 484 * ----------------------------------------------------------------- 485 */ 486 487 SUNDIALS_EXPORT int CVSpilsSetPrecTypeB(void *cvode_mem, int which, int pretypeB); 488 SUNDIALS_EXPORT int CVSpilsSetGSTypeB(void *cvode_mem, int which, int gstypeB); 489 SUNDIALS_EXPORT int CVSpilsSetEpsLinB(void *cvode_mem, int which, realtype eplifacB); 490 SUNDIALS_EXPORT int CVSpilsSetMaxlB(void *cvode_mem, int which, int maxlB); 491 492 SUNDIALS_EXPORT int CVSpilsSetPreconditionerB(void *cvode_mem, int which, 493 CVSpilsPrecSetupFnB psetB, 494 CVSpilsPrecSolveFnB psolveB); 495 SUNDIALS_EXPORT int CVSpilsSetPreconditionerBS(void *cvode_mem, int which, 496 CVSpilsPrecSetupFnBS psetBS, 497 CVSpilsPrecSolveFnBS psolveBS); 498 499 SUNDIALS_EXPORT int CVSpilsSetJacTimesVecFnB(void *cvode_mem, int which, 500 CVSpilsJacTimesVecFnB jtvB); 501 SUNDIALS_EXPORT int CVSpilsSetJacTimesVecFnBS(void *cvode_mem, int which, 502 CVSpilsJacTimesVecFnBS jtvBS); 503 504 505 #ifdef __cplusplus 506 } 507 #endif 508 509 #endif 510