1 /* 2 * ----------------------------------------------------------------- 3 * Programmer(s): Allan Taylor, Alan Hindmarsh, Radu Serban, and 4 * Aaron Collier @ LLNL 5 * ----------------------------------------------------------------- 6 * SUNDIALS Copyright Start 7 * Copyright (c) 2002-2020, Lawrence Livermore National Security 8 * and Southern Methodist University. 9 * All rights reserved. 10 * 11 * See the top-level LICENSE and NOTICE files for details. 12 * 13 * SPDX-License-Identifier: BSD-3-Clause 14 * SUNDIALS Copyright End 15 * ----------------------------------------------------------------- 16 * KINSOL solver module header file (private version) 17 * ----------------------------------------------------------------- 18 */ 19 20 #ifndef _KINSOL_IMPL_H 21 #define _KINSOL_IMPL_H 22 23 #include <stdarg.h> 24 25 #include <kinsol/kinsol.h> 26 27 #ifdef __cplusplus /* wrapper to enable C++ usage */ 28 extern "C" { 29 #endif 30 31 /* 32 * ================================================================= 33 * M A I N S O L V E R M E M O R Y B L O C K 34 * ================================================================= 35 */ 36 37 /* KINSOL default constants */ 38 39 #define PRINTFL_DEFAULT 0 40 #define MXITER_DEFAULT 200 41 #define MXNBCF_DEFAULT 10 42 #define MSBSET_DEFAULT 10 43 #define MSBSET_SUB_DEFAULT 5 44 45 #define OMEGA_MIN RCONST(0.00001) 46 #define OMEGA_MAX RCONST(0.9) 47 48 /* 49 * ----------------------------------------------------------------- 50 * Types : struct KINMemRec and struct *KINMem 51 * ----------------------------------------------------------------- 52 * A variable declaration of type struct *KINMem denotes a 53 * pointer to a data structure of type struct KINMemRec. The 54 * KINMemRec structure contains numerous fields that must be 55 * accessible by KINSOL solver module routines. 56 * ----------------------------------------------------------------- 57 */ 58 59 typedef struct KINMemRec { 60 61 realtype kin_uround; /* machine epsilon (or unit roundoff error) 62 (defined in sundials_types.h) */ 63 64 /* problem specification data */ 65 66 KINSysFn kin_func; /* nonlinear system function implementation */ 67 void *kin_user_data; /* work space available to func routine */ 68 realtype kin_fnormtol; /* stopping tolerance on L2-norm of function 69 value */ 70 realtype kin_scsteptol; /* scaled step length tolerance */ 71 int kin_globalstrategy; /* choices are KIN_NONE, KIN_LINESEARCH 72 KIN_PICARD and KIN_FP */ 73 int kin_printfl; /* level of verbosity of output */ 74 long int kin_mxiter; /* maximum number of nonlinear iterations */ 75 long int kin_msbset; /* maximum number of nonlinear iterations that 76 may be performed between calls to the 77 linear solver setup routine (lsetup) */ 78 long int kin_msbset_sub; /* subinterval length for residual monitoring */ 79 long int kin_mxnbcf; /* maximum number of beta condition failures */ 80 int kin_etaflag; /* choices are KIN_ETACONSTANT, KIN_ETACHOICE1 81 and KIN_ETACHOICE2 */ 82 booleantype kin_noMinEps; /* flag controlling whether or not the value 83 of eps is bounded below */ 84 booleantype kin_constraintsSet; /* flag indicating if constraints are being 85 used */ 86 booleantype kin_jacCurrent; /* flag indicating if the Jacobian info. 87 used by the linear solver is current */ 88 booleantype kin_callForcingTerm; /* flag set if using either KIN_ETACHOICE1 89 or KIN_ETACHOICE2 */ 90 booleantype kin_noResMon; /* flag indicating if the nonlinear 91 residual monitoring scheme should be 92 used */ 93 booleantype kin_retry_nni; /* flag indicating if nonlinear iteration 94 should be retried (set by residual 95 monitoring algorithm) */ 96 booleantype kin_update_fnorm_sub; /* flag indicating if the fnorm associated 97 with the subinterval needs to be 98 updated (set by residual monitoring 99 algorithm) */ 100 101 realtype kin_mxnewtstep; /* maximum allowable scaled step length */ 102 realtype kin_mxnstepin; /* input (or preset) value for mxnewtstep */ 103 realtype kin_sqrt_relfunc; /* relative error bound for func(u) */ 104 realtype kin_stepl; /* scaled length of current step */ 105 realtype kin_stepmul; /* step scaling factor */ 106 realtype kin_eps; /* current value of eps */ 107 realtype kin_eta; /* current value of eta */ 108 realtype kin_eta_gamma; /* gamma value used in eta calculation 109 (choice #2) */ 110 realtype kin_eta_alpha; /* alpha value used in eta calculation 111 (choice #2) */ 112 booleantype kin_noInitSetup; /* flag controlling whether or not the KINSol 113 routine makes an initial call to the 114 linear solver setup routine (lsetup) */ 115 realtype kin_sthrsh; /* threshold value for calling the linear 116 solver setup routine */ 117 118 /* counters */ 119 120 long int kin_nni; /* number of nonlinear iterations */ 121 long int kin_nfe; /* number of calls made to func routine */ 122 long int kin_nnilset; /* value of nni counter when the linear solver 123 setup was last called */ 124 long int kin_nnilset_sub; /* value of nni counter when the linear solver 125 setup was last called (subinterval) */ 126 long int kin_nbcf; /* number of times the beta-condition could not 127 be met in KINLineSearch */ 128 long int kin_nbktrk; /* number of backtracks performed by 129 KINLineSearch */ 130 long int kin_ncscmx; /* number of consecutive steps of size 131 mxnewtstep taken */ 132 133 /* vectors */ 134 135 N_Vector kin_uu; /* solution vector/current iterate (initially 136 contains initial guess, but holds approximate 137 solution upon completion if no errors occurred) */ 138 N_Vector kin_unew; /* next iterate (unew = uu+pp) */ 139 N_Vector kin_fval; /* vector containing result of nonlinear system 140 function evaluated at a given iterate 141 (fval = func(uu)) */ 142 N_Vector kin_gval; /* vector containing result of the fixed point 143 function evaluated at a given iterate; 144 used in KIN_PICARD strategy only. 145 (gval = uu - L^{-1}fval(uu)) */ 146 N_Vector kin_uscale; /* iterate scaling vector */ 147 N_Vector kin_fscale; /* fval scaling vector */ 148 N_Vector kin_pp; /* incremental change vector (pp = unew-uu) */ 149 N_Vector kin_constraints; /* constraints vector */ 150 N_Vector kin_vtemp1; /* scratch vector #1 */ 151 N_Vector kin_vtemp2; /* scratch vector #2 */ 152 153 /* space requirements for AA, Broyden and NLEN */ 154 N_Vector kin_fold_aa; /* vector needed for AA, Broyden, and NLEN */ 155 N_Vector kin_gold_aa; /* vector needed for AA, Broyden, and NLEN */ 156 N_Vector *kin_df_aa; /* vector array needed for AA, Broyden, and NLEN */ 157 N_Vector *kin_dg_aa; /* vector array needed for AA, Broyden and NLEN */ 158 N_Vector *kin_q_aa; /* vector array needed for AA */ 159 realtype kin_beta_aa; /* beta damping parameter for AA */ 160 realtype *kin_gamma_aa; /* array of size maa used in AA */ 161 realtype *kin_R_aa; /* array of size maa*maa used in AA */ 162 long int *kin_ipt_map; /* array of size maa used in AA */ 163 long int kin_m_aa; /* parameter for AA, Broyden or NLEN */ 164 booleantype kin_aamem_aa; /* sets additional memory needed for Anderson Acc */ 165 booleantype kin_setstop_aa; /* determines whether user will set stopping criterion */ 166 booleantype kin_damping_aa; /* flag to apply damping in AA */ 167 realtype *kin_cv; /* scalar array for fused vector operations */ 168 N_Vector *kin_Xv; /* vector array for fused vector operations */ 169 170 /* space requirements for vector storage */ 171 172 sunindextype kin_lrw1; /* number of realtype-sized memory blocks needed 173 for a single N_Vector */ 174 sunindextype kin_liw1; /* number of int-sized memory blocks needed for 175 a single N_Vecotr */ 176 long int kin_lrw; /* total number of realtype-sized memory blocks 177 needed for all KINSOL work vectors */ 178 long int kin_liw; /* total number of int-sized memory blocks needed 179 for all KINSOL work vectors */ 180 181 /* linear solver data */ 182 183 /* function prototypes (pointers) */ 184 185 int (*kin_linit)(struct KINMemRec *kin_mem); 186 187 int (*kin_lsetup)(struct KINMemRec *kin_mem); 188 189 int (*kin_lsolve)(struct KINMemRec *kin_mem, N_Vector xx, N_Vector bb, 190 realtype *sJpnorm, realtype *sFdotJp); 191 192 int (*kin_lfree)(struct KINMemRec *kin_mem); 193 194 booleantype kin_inexact_ls; /* flag set by the linear solver module 195 (in linit) indicating whether this is an 196 iterative linear solver (SUNTRUE), or a direct 197 linear solver (SUNFALSE) */ 198 199 void *kin_lmem; /* pointer to linear solver memory block */ 200 201 realtype kin_fnorm; /* value of L2-norm of fscale*fval */ 202 realtype kin_f1norm; /* f1norm = 0.5*(fnorm)^2 */ 203 realtype kin_sFdotJp; /* value of scaled F(u) vector (fscale*fval) 204 dotted with scaled J(u)*pp vector (set by lsolve) */ 205 realtype kin_sJpnorm; /* value of L2-norm of fscale*(J(u)*pp) 206 (set by lsolve) */ 207 208 realtype kin_fnorm_sub; /* value of L2-norm of fscale*fval (subinterval) */ 209 booleantype kin_eval_omega; /* flag indicating that omega must be evaluated. */ 210 realtype kin_omega; /* constant value for real scalar used in test to 211 determine if reduction of norm of nonlinear 212 residual is sufficient. Unless a valid constant 213 value is specified by the user, omega is estimated 214 from omega_min and omega_max at each iteration. */ 215 realtype kin_omega_min; /* lower bound on omega */ 216 realtype kin_omega_max; /* upper bound on omega */ 217 218 /* 219 * ----------------------------------------------------------------- 220 * Note: The KINLineSearch subroutine scales the values of the 221 * variables sFdotJp and sJpnorm by a factor rl (lambda) that is 222 * chosen by the line search algorithm such that the sclaed Newton 223 * step satisfies the following conditions: 224 * 225 * F(u_k+1) <= F(u_k) + alpha*(F(u_k)^T * J(u_k))*p*rl 226 * 227 * F(u_k+1) >= F(u_k) + beta*(F(u_k)^T * J(u_k))*p*rl 228 * 229 * where alpha = 1.0e-4, beta = 0.9, u_k+1 = u_k + rl*p, 230 * 0 < rl <= 1, J denotes the system Jacobian, and F represents 231 * the nonliner system function. 232 * ----------------------------------------------------------------- 233 */ 234 235 booleantype kin_MallocDone; /* flag indicating if KINMalloc has been 236 called yet */ 237 238 /* message files */ 239 /*------------------------------------------- 240 Error handler function and error ouput file 241 -------------------------------------------*/ 242 243 KINErrHandlerFn kin_ehfun; /* Error messages are handled by ehfun */ 244 void *kin_eh_data; /* dats pointer passed to ehfun */ 245 FILE *kin_errfp; /* KINSOL error messages are sent to errfp */ 246 247 KINInfoHandlerFn kin_ihfun; /* Info messages are handled by ihfun */ 248 void *kin_ih_data; /* dats pointer passed to ihfun */ 249 FILE *kin_infofp; /* where KINSol info messages are sent */ 250 251 } *KINMem; 252 253 /* 254 * ================================================================= 255 * I N T E R F A C E T O L I N E A R S O L V E R 256 * ================================================================= 257 */ 258 259 /* 260 * ----------------------------------------------------------------- 261 * Function : int (*kin_linit)(KINMem kin_mem) 262 * ----------------------------------------------------------------- 263 * kin_linit initializes solver-specific data structures (including 264 * variables used as counters or for storing statistical information), 265 * but system memory allocation should be done by the subroutine 266 * that actually initializes the environment for liner solver 267 * package. If the linear system is to be preconditioned, then the 268 * variable setupNonNull (type booleantype) should be set to SUNTRUE 269 * (predefined constant) and the kin_lsetup routine should be 270 * appropriately defined. 271 * 272 * kinmem pointer to an internal memory block allocated during 273 * prior calls to KINCreate and KINMalloc 274 * 275 * If the necessary variables have been successfully initialized, 276 * then the kin_linit function should return 0 (zero). Otherwise, 277 * the subroutine should indicate a failure has occurred by 278 * returning a non-zero integer value. 279 * ----------------------------------------------------------------- 280 */ 281 282 /* 283 * ----------------------------------------------------------------- 284 * Function : int (*kin_lsetup)(KINMem kin_mem) 285 * ----------------------------------------------------------------- 286 * kin_lsetup interfaces with the user-supplied pset subroutine (the 287 * preconditioner setup routine), and updates relevant variable 288 * values (see KINSpgmrSetup/KINSpbcgSetup). Simply stated, the 289 * kin_lsetup routine prepares the linear solver for a subsequent 290 * call to the user-supplied kin_lsolve function. 291 * 292 * kinmem pointer to an internal memory block allocated during 293 * prior calls to KINCreate and KINMalloc 294 * 295 * If successful, the kin_lsetup routine should return 0 (zero). 296 * Otherwise it should return a non-zero value. 297 * ----------------------------------------------------------------- 298 */ 299 300 /* 301 * ----------------------------------------------------------------- 302 * Function : int (*kin_lsolve)(KINMem kin_mem, N_Vector xx, 303 * N_Vector bb, realtype *sJpnorm, realtype *sFdotJp) 304 * ----------------------------------------------------------------- 305 * kin_lsolve interfaces with the subroutine implementing the 306 * numerical method to be used to solve the linear system J*xx = bb, 307 * and must increment the relevant counter variable values in 308 * addition to computing certain values used by the global strategy 309 * and forcing term routines (see KINInexactNewton, KINLineSearch, 310 * KINForcingTerm, and KINSpgmrSolve/KINSpbcgSolve). 311 * 312 * kinmem pointer to an internal memory block allocated during 313 * prior calls to KINCreate and KINMalloc 314 * 315 * xx vector (type N_Vector) set to initial guess by kin_lsolve 316 * routine prior to calling the linear solver, but which upon 317 * return contains an approximate solution of the linear 318 * system J*xx = bb, where J denotes the system Jacobian 319 * 320 * bb vector (type N_Vector) set to -func(u) (negative of the 321 * value of the system function evaluated at the current 322 * iterate) by KINLinSolDrv before kin_lsolve is called 323 * 324 * sJpnorm holds the value of the L2-norm (Euclidean norm) of 325 * fscale*(J(u)*pp) upon return 326 * 327 * sFdotJp holds the value of the scaled F(u) (fscale*F) dotted 328 * with the scaled J(u)*pp vector upon return 329 * 330 * If successful, the kin_lsolve routine should return 0 (zero). 331 * Otherwise it should return a positive value if a re-evaluation 332 * of the lsetup function could recover, or a negative value if 333 * no such recovery is possible. 334 * ----------------------------------------------------------------- 335 */ 336 337 /* 338 * ----------------------------------------------------------------- 339 * Function : int (*kin_lfree)(KINMem kin_mem) 340 * ----------------------------------------------------------------- 341 * kin_lfree is called by KINFree and should free (deallocate) all 342 * system memory resources allocated for the linear solver module 343 * (see KINSpgmrFree/KINSpbcgFree). It should return 0 upon 344 * success, nonzero on failure. 345 * 346 * kinmem pointer to an internal memory block allocated during 347 * prior calls to KINCreate and KINMalloc 348 * ----------------------------------------------------------------- 349 */ 350 351 /* 352 * ================================================================= 353 * K I N S O L I N T E R N A L F U N C T I O N S 354 * ================================================================= 355 */ 356 357 358 /* High level error handler */ 359 360 void KINProcessError(KINMem kin_mem, 361 int error_code, const char *module, const char *fname, 362 const char *msgfmt, ...); 363 364 /* Prototype of internal errHandler function */ 365 366 void KINErrHandler(int error_code, const char *module, const char *function, 367 char *msg, void *user_data); 368 369 370 /* High level info handler */ 371 372 void KINPrintInfo(KINMem kin_mem, 373 int info_code, const char *module, const char *fname, 374 const char *msgfmt, ...); 375 376 /* Prototype of internal infoHandler function */ 377 378 void KINInfoHandler(const char *module, const char *function, 379 char *msg, void *user_data); 380 381 /* 382 * ================================================================= 383 * K I N S O L E R R O R M E S S A G E S 384 * ================================================================= 385 */ 386 387 #define MSG_MEM_FAIL "A memory request failed." 388 #define MSG_NO_MEM "kinsol_mem = NULL illegal." 389 #define MSG_BAD_NVECTOR "A required vector operation is not implemented." 390 #define MSG_FUNC_NULL "func = NULL illegal." 391 #define MSG_NO_MALLOC "Attempt to call before KINMalloc illegal." 392 393 #define MSG_BAD_PRINTFL "Illegal value for printfl." 394 #define MSG_BAD_MXITER "Illegal value for mxiter." 395 #define MSG_BAD_MSBSET "Illegal msbset < 0." 396 #define MSG_BAD_MSBSETSUB "Illegal msbsetsub < 0." 397 #define MSG_BAD_ETACHOICE "Illegal value for etachoice." 398 #define MSG_BAD_ETACONST "eta out of range." 399 #define MSG_BAD_GAMMA "gamma out of range." 400 #define MSG_BAD_ALPHA "alpha out of range." 401 #define MSG_BAD_MXNEWTSTEP "Illegal mxnewtstep < 0." 402 #define MSG_BAD_RELFUNC "relfunc < 0 illegal." 403 #define MSG_BAD_FNORMTOL "fnormtol < 0 illegal." 404 #define MSG_BAD_SCSTEPTOL "scsteptol < 0 illegal." 405 #define MSG_BAD_MXNBCF "mxbcf < 0 illegal." 406 #define MSG_BAD_CONSTRAINTS "Illegal values in constraints vector." 407 #define MSG_BAD_OMEGA "scalars < 0 illegal." 408 #define MSG_BAD_MAA "maa < 0 illegal." 409 #define MSG_ZERO_MAA "maa = 0 illegal." 410 411 #define MSG_LSOLV_NO_MEM "The linear solver memory pointer is NULL." 412 #define MSG_UU_NULL "uu = NULL illegal." 413 #define MSG_BAD_GLSTRAT "Illegal value for global strategy." 414 #define MSG_BAD_USCALE "uscale = NULL illegal." 415 #define MSG_USCALE_NONPOSITIVE "uscale has nonpositive elements." 416 #define MSG_BAD_FSCALE "fscale = NULL illegal." 417 #define MSG_FSCALE_NONPOSITIVE "fscale has nonpositive elements." 418 #define MSG_CONSTRAINTS_NOTOK "Constraints not allowed with fixed point or Picard iterations" 419 #define MSG_INITIAL_CNSTRNT "Initial guess does NOT meet constraints." 420 #define MSG_LINIT_FAIL "The linear solver's init routine failed." 421 422 #define MSG_SYSFUNC_FAILED "The system function failed in an unrecoverable manner." 423 #define MSG_SYSFUNC_FIRST "The system function failed at the first call." 424 #define MSG_LSETUP_FAILED "The linear solver's setup function failed in an unrecoverable manner." 425 #define MSG_LSOLVE_FAILED "The linear solver's solve function failed in an unrecoverable manner." 426 #define MSG_LINSOLV_NO_RECOVERY "The linear solver's solve function failed recoverably, but the Jacobian data is already current." 427 #define MSG_LINESEARCH_NONCONV "The line search algorithm was unable to find an iterate sufficiently distinct from the current iterate." 428 #define MSG_LINESEARCH_BCFAIL "The line search algorithm was unable to satisfy the beta-condition for nbcfails iterations." 429 #define MSG_MAXITER_REACHED "The maximum number of iterations was reached before convergence." 430 #define MSG_MXNEWT_5X_EXCEEDED "Five consecutive steps have been taken that satisfy a scaled step length test." 431 #define MSG_SYSFUNC_REPTD "Unable to correct repeated recoverable system function errors." 432 #define MSG_NOL_FAIL "Unable to find user's Linear Jacobian, which is required for the KIN_PICARD Strategy" 433 434 /* 435 * ================================================================= 436 * K I N S O L I N F O M E S S A G E S 437 * ================================================================= 438 */ 439 440 #define INFO_RETVAL "Return value: %d" 441 #define INFO_ADJ "no. of lambda adjustments = %ld" 442 443 #if defined(SUNDIALS_EXTENDED_PRECISION) 444 445 #define INFO_NNI "nni = %4ld nfe = %6ld fnorm = %26.16Lg" 446 #define INFO_TOL "scsteptol = %12.3Lg fnormtol = %12.3Lg" 447 #define INFO_FMAX "scaled f norm (for stopping) = %12.3Lg" 448 #define INFO_PNORM "pnorm = %12.4Le" 449 #define INFO_PNORM1 "(ivio=1) pnorm = %12.4Le" 450 #define INFO_FNORM "fnorm(L2) = %20.8Le" 451 #define INFO_LAM "min_lam = %11.4Le f1norm = %11.4Le pnorm = %11.4Le" 452 #define INFO_ALPHA "fnorm = %15.8Le f1norm = %15.8Le alpha_cond = %15.8Le lam = %15.8Le" 453 #define INFO_BETA "f1norm = %15.8Le beta_cond = %15.8Le lam = %15.8Le" 454 #define INFO_ALPHABETA "f1norm = %15.8Le alpha_cond = %15.8Le beta_cond = %15.8Le lam = %15.8Le" 455 456 #elif defined(SUNDIALS_DOUBLE_PRECISION) 457 458 #define INFO_NNI "nni = %4ld nfe = %6ld fnorm = %26.16lg" 459 #define INFO_TOL "scsteptol = %12.3lg fnormtol = %12.3lg" 460 #define INFO_FMAX "scaled f norm (for stopping) = %12.3lg" 461 #define INFO_PNORM "pnorm = %12.4le" 462 #define INFO_PNORM1 "(ivio=1) pnorm = %12.4le" 463 #define INFO_FNORM "fnorm(L2) = %20.8le" 464 #define INFO_LAM "min_lam = %11.4le f1norm = %11.4le pnorm = %11.4le" 465 #define INFO_ALPHA "fnorm = %15.8le f1norm = %15.8le alpha_cond = %15.8le lam = %15.8le" 466 #define INFO_BETA "f1norm = %15.8le beta_cond = %15.8le lam = %15.8le" 467 #define INFO_ALPHABETA "f1norm = %15.8le alpha_cond = %15.8le beta_cond = %15.8le lam = %15.8le" 468 469 #else 470 471 #define INFO_NNI "nni = %4ld nfe = %6ld fnorm = %26.16g" 472 #define INFO_TOL "scsteptol = %12.3g fnormtol = %12.3g" 473 #define INFO_FMAX "scaled f norm (for stopping) = %12.3g" 474 #define INFO_PNORM "pnorm = %12.4e" 475 #define INFO_PNORM1 "(ivio=1) pnorm = %12.4e" 476 #define INFO_FNORM "fnorm(L2) = %20.8e" 477 #define INFO_LAM "min_lam = %11.4e f1norm = %11.4e pnorm = %11.4e" 478 #define INFO_ALPHA "fnorm = %15.8e f1norm = %15.8e alpha_cond = %15.8e lam = %15.8e" 479 #define INFO_BETA "f1norm = %15.8e beta_cond = %15.8e lam = %15.8e" 480 #define INFO_ALPHABETA "f1norm = %15.8e alpha_cond = %15.8e beta_cond = %15.8e lam = %15.8e" 481 482 #endif 483 484 485 #ifdef __cplusplus 486 } 487 #endif 488 489 #endif 490