1 /* 2 * ----------------------------------------------------------------- 3 * $Revision: 1.5 $ 4 * $Date: 2007/10/26 21:51:29 $ 5 * ----------------------------------------------------------------- 6 * Programmer: Radu Serban @ LLNL 7 * ----------------------------------------------------------------- 8 * Copyright (c) 2006, The Regents of the University of California. 9 * Produced at the Lawrence Livermore National Laboratory. 10 * All rights reserved. 11 * For details, see the LICENSE file. 12 * ----------------------------------------------------------------- 13 * Implementation header file for the main CPODES integrator. 14 * ----------------------------------------------------------------- 15 */ 16 17 #ifndef _CPODES_IMPL_H 18 #define _CPODES_IMPL_H 19 20 #ifdef __cplusplus /* wrapper to enable C++ usage */ 21 extern "C" { 22 #endif 23 24 #include <cpodes/cpodes.h> 25 26 /* 27 * ================================================================= 28 * M A I N I N T E G R A T O R M E M O R Y B L O C K 29 * ================================================================= 30 */ 31 32 /* Basic LMM constants */ 33 34 #define ADAMS_Q_MAX 12 /* max value of q for ADAMS */ 35 #define BDF_Q_MAX 5 /* max value of q for BDF */ 36 #define Q_MAX ADAMS_Q_MAX /* max value of q for either LMM */ 37 #define L_MAX (Q_MAX+1) /* max value of L for either LMM */ 38 #define NUM_TESTS 5 /* number of error test quantities */ 39 40 /* 41 * ----------------------------------------------------------------- 42 * Types : struct CPodeMemRec, CPodeMem 43 * ----------------------------------------------------------------- 44 * The type CPodeMem is type pointer to struct CPodeMemRec. 45 * This structure contains fields to keep track of problem state. 46 * ----------------------------------------------------------------- 47 */ 48 49 typedef struct CPodeMemRec { 50 51 realtype cp_uround; /* machine unit roundoff */ 52 53 /*-------------------------- 54 Problem Specification Data 55 --------------------------*/ 56 57 int cp_ode_type; /* ODE type: ode = CP_IMPL or CP_EXPL */ 58 CPResFn cp_fi; /* F(t,y'(t),y(t)) = 0 */ 59 CPRhsFn cp_fe; /* y' = f(t,y(t)) */ 60 void *cp_f_data; /* user pointer passed to f */ 61 62 int cp_lmm_type; /* lmm_type = CP_ADAMS or CP_BDF */ 63 int cp_nls_type; /* nls_type = CP_FUNCTIONAL or CP_NEWTON */ 64 int cp_tol_type; /* tol_type = CP_SS or CP_SV or CP_WF */ 65 66 realtype cp_reltol; /* relative tolerance */ 67 realtype cp_Sabstol; /* scalar absolute tolerance */ 68 N_Vector cp_Vabstol; /* vector absolute tolerance */ 69 CPEwtFn cp_efun; /* function to set ewt */ 70 void *cp_e_data; /* user pointer passed to efun */ 71 72 /*----------------------- 73 Nordsieck History Array 74 -----------------------*/ 75 76 /* 77 * Nordsieck array, of size N x (q+1). 78 * zn[j] is a vector of length N (j=0,...,q) 79 * zn[j] = [1/factorial(j)] * h^j * (jth derivative of the interpolating polynomial) 80 */ 81 82 N_Vector cp_zn[L_MAX]; 83 84 /*------------------------------------------------ 85 Other vectors of same length as the state vector 86 ------------------------------------------------*/ 87 88 N_Vector cp_ewt; /* error weight vector. */ 89 N_Vector cp_y; /* work space for y vector (= yout) */ 90 N_Vector cp_yp; /* work space for yp vector (= ypout) */ 91 N_Vector cp_acor; /* accumulated correction acor = y_n(m)-y_n(0) */ 92 N_Vector cp_tempv; /* temporary storage vector */ 93 N_Vector cp_ftemp; /* temporary storage vector */ 94 95 /*--------- 96 Step Data 97 ---------*/ 98 99 int cp_q; /* current order */ 100 int cp_qprime; /* order to be used on the next step */ 101 int cp_next_q; /* order to be used on the next step */ 102 int cp_qwait; /* no. of steps to wait before a change in q */ 103 int cp_L; /* L = q + 1 */ 104 105 realtype cp_hin; /* initial step size */ 106 realtype cp_h; /* current step size */ 107 realtype cp_hprime; /* step size to be used on the next step */ 108 realtype cp_next_h; /* step size to be used on the next step */ 109 realtype cp_eta; /* eta = hprime / h */ 110 realtype cp_hscale; /* value of h used in zn */ 111 realtype cp_tn; /* current internal value of t */ 112 realtype cp_tretlast; /* value of tret last returned by CPode */ 113 114 realtype cp_tau[L_MAX+1]; /* array of previous q+1 successful step sizes */ 115 realtype cp_tq[NUM_TESTS+1]; /* array of test quantities */ 116 realtype cp_l[L_MAX]; /* coefficients of Lambda(x) (degree q poly) */ 117 realtype cp_p[L_MAX]; /* coefficients of Phi(x) (degree q poly) */ 118 119 realtype cp_rl1; /* the scalar 1/l[1] */ 120 realtype cp_gamma; /* gamma = h * rl1 */ 121 realtype cp_gammap; /* gamma at the last setup call */ 122 realtype cp_gamrat; /* gamma / gammap */ 123 124 realtype cp_crate; /* estimated corrector convergence rate */ 125 realtype cp_acnrm; /* ||acor||_wrms */ 126 realtype cp_nlscoef; /* coeficient in nonlinear convergence test */ 127 int cp_mnewt; /* Newton iteration counter */ 128 129 /*--------------- 130 Projection Data 131 ---------------*/ 132 133 int cp_proj_type; /* user vs. internal projection algorithm */ 134 booleantype cp_proj_enabled; /* is projection enabled? */ 135 booleantype cp_applyProj; /* was projection performed at current step? */ 136 137 long int cp_proj_freq; /* projection frequency */ 138 long int cp_nstlprj; /* step number of last projection */ 139 140 int cp_proj_norm; /* type of projection norm (L2 or WRMS) */ 141 int cp_cnstr_type; /* type of constraints (lin. or nonlin.) */ 142 CPCnstrFn cp_cfun; /* 0 = c(t,y(t)) */ 143 void *cp_c_data; /* user pointer passed to cfun */ 144 145 CPProjFn cp_pfun; /* function to perform projection */ 146 void *cp_p_data; /* user pointer passed to pfun */ 147 148 realtype cp_prjcoef; /* coefficient in projection convergence test */ 149 150 N_Vector cp_acorP; /* projection correction (length N) */ 151 N_Vector cp_errP; /* projected error estimate (length N) */ 152 N_Vector cp_yC; /* saved corrected state (length N) */ 153 154 long int cp_lsetupP_freq; /* frequency of cnstr. Jacobian evaluation */ 155 long int cp_nstlsetP; /* step number of last lsetupP call */ 156 157 realtype cp_crateP; /* estimated conv. rate in nonlin. projection */ 158 int cp_maxcorP; /* maximum number of nonlin. proj. iterations */ 159 160 booleantype cp_project_err; /* should we project the error estimate? */ 161 162 booleantype cp_test_cnstr; /* test constraint norm before projection? */ 163 N_Vector cp_ctol; /* vector of constraint "absolute tolerances" */ 164 165 N_Vector cp_ctemp; /* temporary vectors (length M) */ 166 N_Vector cp_tempvP1; 167 N_Vector cp_tempvP2; 168 169 booleantype cp_first_proj; /* is this the first time we project? */ 170 171 long int cp_nproj; /* number of projection steps performed */ 172 long int cp_nprf; /* number of projection failures */ 173 long int cp_nce; /* number of calls to cfun */ 174 long int cp_nsetupsP; /* number of calls to lsetupP */ 175 int cp_maxnpf; /* maximum number of projection failures */ 176 177 /*----------------------- 178 Quadrature Related Data 179 -----------------------*/ 180 181 booleantype cp_quadr; /* are we integrating quadratures? */ 182 183 CPQuadFn cp_qfun; /* function defining the integrand */ 184 void *cp_q_data; /* user pointer passed to fQ */ 185 int cp_tol_typeQ; /* tol_typeQ = CP_SS or CP_SV */ 186 booleantype cp_errconQ; /* are quadratures included in error test? */ 187 188 realtype cp_reltolQ; /* relative tolerance for quadratures */ 189 realtype cp_SabstolQ; /* scalar absolute tolerance for quadratures */ 190 N_Vector cp_VabstolQ; /* vector absolute tolerance for quadratures */ 191 192 N_Vector cp_znQ[L_MAX]; /* Nordsieck arrays for sensitivities */ 193 N_Vector cp_ewtQ; /* error weight vector for quadratures */ 194 N_Vector cp_yQ; /* Unlike y, yQ is not allocated by the user */ 195 N_Vector cp_acorQ; /* acorQ = yQ_n(m) - yQ_n(0) */ 196 N_Vector cp_tempvQ; /* temporary storage vector (~ tempv) */ 197 198 realtype cp_acnrmQ; /* acnrmQ = ||acorQ||_WRMS */ 199 200 long int cp_lrw1Q; /* no. of realtype words in 1 N_Vector yQ */ 201 long int cp_liw1Q; /* no. of integer words in 1 N_Vector yQ */ 202 203 int cp_qmax_allocQ; /* qmax used when allocating quad. memory */ 204 205 booleantype cp_VabstolQMallocDone; /* did we allocate memory for abstolQ? */ 206 booleantype cp_quadMallocDone; /* was quadrature memory allocated? */ 207 208 long int cp_nqe; /* number of calls to qfun */ 209 long int cp_netfQ; /* number of quadr. error test failures */ 210 211 /*----------------- 212 Tstop information 213 -----------------*/ 214 215 booleantype cp_tstopset; 216 booleantype cp_istop; 217 realtype cp_tstop; 218 219 /*------ 220 Limits 221 ------*/ 222 223 int cp_qmax; /* q <= qmax */ 224 long int cp_mxstep; /* maximum number of internal steps for one user call */ 225 int cp_maxcor; /* maximum number of corrector iterations for the */ 226 /* solution of the nonlinear equation */ 227 int cp_mxhnil; /* maximum number of warning messages issued to the */ 228 /* user that t + h == t for the next internal step */ 229 int cp_maxnef; /* maximum number of error test failures */ 230 int cp_maxncf; /* maximum number of nonlinear convergence failures */ 231 232 realtype cp_hmin; /* |h| >= hmin */ 233 realtype cp_hmax_inv; /* |h| <= 1/hmax_inv */ 234 realtype cp_etamax; /* eta <= etamax */ 235 236 /*-------- 237 Counters 238 --------*/ 239 240 long int cp_nst; /* number of internal steps taken */ 241 long int cp_nfe; /* number of f calls */ 242 long int cp_ncfn; /* number of corrector convergence failures */ 243 long int cp_netf; /* number of error test failures */ 244 long int cp_nni; /* number of nonlinear iterations performed */ 245 long int cp_nsetups; /* number of calls to lsetup */ 246 int cp_nhnil; /* number of messages saying that t+h==t */ 247 248 realtype cp_etaqm1; /* ratio of new to old h for order q-1 */ 249 realtype cp_etaq; /* ratio of new to old h for order q */ 250 realtype cp_etaqp1; /* ratio of new to old h for order q+1 */ 251 252 /*---------------------------- 253 Space requirements for CPODES 254 ----------------------------*/ 255 256 long int cp_lrw1; /* no. of realtype words in 1 state N_Vector */ 257 long int cp_liw1; /* no. of integer words in 1 state N_Vector */ 258 long int cp_lrw2; /* no. of realtype words in 1 cnstr. N_Vector */ 259 long int cp_liw2; /* no. of integer words in 1 cnstr. N_Vector */ 260 long int cp_lrw; /* no. of realtype words in CPODES work vectors*/ 261 long int cp_liw; /* no. of integer words in CPODES work vectors */ 262 263 /*--------------------------------------- 264 Implicit Integration Linear Solver Data 265 ---------------------------------------*/ 266 267 /* Linear Solver functions to be called */ 268 int (*cp_linit)(struct CPodeMemRec *cp_mem); 269 int (*cp_lsetup)(struct CPodeMemRec *cp_mem, int convfail, 270 N_Vector yP, N_Vector ypP, N_Vector fctP, 271 booleantype *jcurPtr, 272 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 273 int (*cp_lsolve)(struct CPodeMemRec *cp_mem, N_Vector b, N_Vector weight, 274 N_Vector yC, N_Vector ypC, N_Vector fctC); 275 void (*cp_lfree)(struct CPodeMemRec *cp_mem); 276 277 /* Linear Solver specific memory */ 278 void *cp_lmem; 279 280 /* Does lsetup do anything? */ 281 booleantype cp_lsetup_exists; 282 283 /*----------------------------- 284 Projection Linear Solver Data 285 -----------------------------*/ 286 287 /* Linear Solver functions to be called */ 288 int (*cp_linitP)(struct CPodeMemRec *cp_mem); 289 int (*cp_lsetupP)(struct CPodeMemRec *cp_mem, 290 N_Vector y, N_Vector cy, 291 N_Vector c_tmp1, N_Vector c_tmp2, N_Vector s_tmp1); 292 int (*cp_lsolveP)(struct CPodeMemRec *cp_mem, N_Vector b, N_Vector x, 293 N_Vector y, N_Vector cy, 294 N_Vector c_tmp1, N_Vector s_tmp1); 295 void (*cp_lmultP)(struct CPodeMemRec *cp_mem, N_Vector x, N_Vector Gx); 296 void (*cp_lfreeP)(struct CPodeMemRec *cp_mem); 297 298 /* Linear Solver specific memory */ 299 void *cp_lmemP; 300 301 /* Does lsetupP do anything? */ 302 booleantype cp_lsetupP_exists; 303 304 /*------------ 305 Saved Values 306 ------------*/ 307 308 int cp_qu; /* last successful q value used */ 309 long int cp_nstlset; /* step number of last lsetup call */ 310 realtype cp_h0u; /* actual initial stepsize */ 311 realtype cp_hu; /* last successful h value used */ 312 realtype cp_saved_tq5; /* saved value of tq[5] */ 313 booleantype cp_jcur; /* Is Jacobian info current? */ 314 realtype cp_tolsf; /* tolerance scale factor */ 315 int cp_qmax_alloc; /* value of qmax used when allocating memory */ 316 int cp_indx_acor; /* index of the zn vector where acor is saved */ 317 318 booleantype cp_MallocDone; 319 booleantype cp_VabstolMallocDone; 320 booleantype cp_projMallocDone; 321 booleantype cp_rootMallocDone; 322 323 324 /*------------------------------------------- 325 Error handler function and error output file 326 -------------------------------------------*/ 327 328 CPErrHandlerFn cp_ehfun; /* Error messages are handled by ehfun */ 329 void *cp_eh_data; /* user pointer passed to ehfun */ 330 FILE *cp_errfp; /* CPODES error messages are sent to errfp */ 331 332 /*------------------------- 333 Stability Limit Detection 334 -------------------------*/ 335 336 booleantype cp_sldeton; /* Is Stability Limit Detection on? */ 337 realtype cp_ssdat[6][4]; /* scaled data array for STALD */ 338 int cp_nscon; /* counter for STALD method */ 339 long int cp_nor; /* counter for number of order reductions */ 340 341 /*---------------- 342 Rootfinding Data 343 ----------------*/ 344 345 booleantype cp_doRootfinding;/* Is rootfinding enabled? */ 346 CPRootFn cp_gfun; /* Function g for roots sought */ 347 int cp_nrtfn; /* number of components of g */ 348 void *cp_g_data; /* pointer to user data for g */ 349 booleantype *cp_gactive; /* flags for active/inactive g functions */ 350 int *cp_iroots; /* array for root information */ 351 int *cp_rootdir; /* array specifying direction of zero-crossing */ 352 353 realtype cp_tlo; /* nearest endpoint of interval in root search */ 354 realtype cp_thi; /* farthest endpoint of interval in root search*/ 355 realtype cp_trout; /* t value returned by rootfinding routine */ 356 realtype *cp_glo; /* saved array of g values at t = tlo */ 357 realtype *cp_ghi; /* saved array of g values at t = thi */ 358 realtype *cp_grout; /* array of g values at t = trout */ 359 realtype cp_toutc; /* copy of tout (if NORMAL mode) */ 360 int cp_taskc; /* copy of parameter task */ 361 int cp_irfnd; /* flag showing whether last step had a root */ 362 long int cp_nge; /* counter for g evaluations */ 363 364 /*------------------------------ 365 Consistent IC calculation data 366 ------------------------------*/ 367 368 realtype cp_icprj_convcoef; 369 realtype cp_icprj_normtol; 370 int cp_icprj_maxrcvr; 371 int cp_icprj_maxiter; 372 N_Vector cp_yy0; 373 N_Vector cp_yp0; 374 375 } *CPodeMem; 376 377 /* 378 * ================================================================= 379 * I N T E R F A C E T O L I N E A R S O L V E R S 380 * F O R I M P L I C I T I N T E G R A T I O N 381 * ================================================================= 382 */ 383 384 /* 385 * ----------------------------------------------------------------- 386 * int (*cp_linit)(CPodeMem cp_mem); 387 * ----------------------------------------------------------------- 388 * The purpose of cp_linit is to complete initializations for a 389 * specific linear solver, such as counters and statistics. 390 * An LInitFn should return 0 if it has successfully initialized the 391 * CPODES linear solver and a negative value otherwise. 392 * If an error does occur, an appropriate message should be sent to 393 * the error handler function. 394 * ----------------------------------------------------------------- 395 */ 396 397 /* 398 * ----------------------------------------------------------------- 399 * int (*cp_lsetup)(CPodeMem cp_mem, int convfail, 400 * N_Vector yP, N_Vector ypP, N_Vector fctP, 401 * booleantype *jcurPtr, 402 * N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 403 * ----------------------------------------------------------------- 404 * The job of cp_lsetup is to prepare the linear solver for subsequent 405 * calls to cp_lsolve. It may recompute Jacobian-related data is it 406 * deems necessary. Its parameters are as follows: 407 * 408 * cp_mem - problem memory pointer of type CPodeMem. See the 409 * typedef earlier in this file. 410 * 411 * convfail - a flag to indicate any problem that occurred during 412 * the solution of the nonlinear equation on the 413 * current time step for which the linear solver is 414 * being used. This flag can be used to help decide 415 * whether the Jacobian data kept by a CPODES linear 416 * solver needs to be updated or not. 417 * Its possible values are documented below. 418 * 419 * NOTE: convfail is UNDEFINED for implicit-form ODE 420 * 421 * yP - the predicted y vector for the current CPODES 422 * internal step (i.e. zn[0]) 423 * 424 * ypP - the predicted y' vector for the current CPODES 425 * internal step (i.e. h*zn[1] for CP_IMPL) 426 * 427 * NOTE: ypP = NULL for explicit-form ODE 428 * 429 * fctP - f(yP) or F(tn, yP, ypP). 430 * 431 * jcurPtr - a pointer to a boolean to be filled in by cp_lsetup. 432 * The function should set *jcurPtr=TRUE if its Jacobian 433 * data is current after the call and should set 434 * *jcurPtr=FALSE if its Jacobian data is not current. 435 * Note: If cp_lsetup calls for re-evaluation of 436 * Jacobian data (based on convfail and CPODES state 437 * data), it should return *jcurPtr=TRUE always; 438 * otherwise an infinite loop can result. 439 * 440 * NOTE: jcurPtr is IGNORED for implicit-form ODE 441 * 442 * tmp1, tmp2, tmp3 - temporary N_Vectors provided for use by cp_lsetup. 443 * 444 * The cp_lsetup routine should return 0 if successful, a positive 445 * value for a recoverable error, and a negative value for an 446 * unrecoverable error. 447 * ----------------------------------------------------------------- 448 */ 449 450 /* 451 * ----------------------------------------------------------------- 452 * int (*cp_lsolve)(CPodeMem cp_mem, N_Vector b, N_Vector weight, 453 * N_Vector yC, N_Vector ypC, N_Vector fctC); 454 * ----------------------------------------------------------------- 455 * cp_lsolve must solve the linear equation P x = b, where 456 * Explicit-form ODE: 457 * P is some approximation to (I - gamma * df/dy) 458 * and the RHS vector b is input. 459 * Implicit-form ODE: 460 * P is some approximation to (dF/dy' + gamma * dF/dy) 461 * and b is the current residual. 462 * Its arguments are as follows: 463 * 464 * cp_mem - problem memory pointer of type CPodeMem. See the big 465 * typedef earlier in this file. 466 * 467 * b - on input the right-hand side of the linear system 468 * to be solved. On output, the solution. 469 * 470 * weight - te current weights in the WRMS norm. 471 * 472 * yC - the solver's current approximation to y(tn) 473 * 474 * ypC - the solver's current approximation to y'(tn) 475 * 476 * NOTE: ypC = NULL for explicit-form ODE 477 * 478 * fctC - f(tn, yC) or F(tn, yC, ypc). 479 * 480 * The solution is to be returned in the vector b. 481 * The cp_lsolve function should return a positive value for a 482 * recoverable error and a negative value for an unrecoverable error. 483 * Success is indicated by a 0 return value. 484 * ----------------------------------------------------------------- 485 */ 486 487 /* 488 * ----------------------------------------------------------------- 489 * void (*cp_lfree)(CPodeMem cp_mem); 490 * ----------------------------------------------------------------- 491 * cp_lfree should free up any memory allocated by the linear 492 * solver. This routine is called once a problem has been 493 * completed and the linear solver is no longer needed. 494 * ----------------------------------------------------------------- 495 */ 496 497 /* 498 * ----------------------------------------------------------------- 499 * Communication between CPODES and a CPODES Linear Solver 500 * ----------------------------------------------------------------- 501 * convfail is passed as an input to cp_lsetup. When dealing with 502 * an ODE in explicit form, this can be used to test whether the Jacobian 503 * must be reevakluated or whether it is enough to use the up-to-date gamma 504 * value. 505 * 506 * convfail can be one of: 507 * 508 * CP_NO_FAILURES : Either this is the first cp_setup call for this 509 * step, or the local error test failed on the 510 * previous attempt at this step (but the Newton 511 * iteration converged). 512 * 513 * CP_FAIL_BAD_J : This value is passed to cp_lsetup if 514 * 515 * (a) The previous Newton corrector iteration 516 * did not converge and the linear solver's 517 * setup routine indicated that its Jacobian- 518 * related data is not current 519 * or 520 * (b) During the previous Newton corrector 521 * iteration, the linear solver's solve routine 522 * failed in a recoverable manner and the 523 * linear solver's setup routine indicated that 524 * its Jacobian-related data is not current. 525 * 526 * CP_FAIL_OTHER : During the current internal step try, the 527 * previous Newton iteration failed to converge 528 * even though the linear solver was using current 529 * Jacobian-related data. 530 * ----------------------------------------------------------------- 531 */ 532 533 534 /* 535 * ================================================================= 536 * I N T E R F A C E T O L I N E A R S O L V E R S 537 * F O R P R O J E C T I O N 538 * ================================================================= 539 */ 540 541 /* 542 * ----------------------------------------------------------------- 543 * int (*cp_linitP)(CPodeMem cp_mem); 544 * ----------------------------------------------------------------- 545 * The purpose of cp_linitP is to complete initializations for a 546 * specific linear solver, such as counters and statistics. 547 * An linit function should return 0 if it has successfully 548 * initialized the CPODES linear solver and a negative value otherwise. 549 * If an error does occur, an appropriate message should be sent to 550 * the error handler function. 551 * ----------------------------------------------------------------- 552 */ 553 554 /* 555 * ----------------------------------------------------------------- 556 * int (*cp_lsetupP)(CPodeMem cp_mem, 557 * N_Vector y, N_Vector cy, 558 * N_Vector c_tmp1, N_Vector c_tmp2, N_Vector s_tmp1); 559 * ----------------------------------------------------------------- 560 * The job of cp_lsetupP is to prepare the linear solver for subsequent 561 * calls to cp_lsolveP. It may recompute Jacobian-related data if it 562 * deems necessary. Its parameters are as follows: 563 * 564 * cp_mem - problem memory pointer of type CPodeMem. See the 565 * typedef earlier in this file. 566 * 567 * y - the corrected y vector for the current CPODES 568 * internal step (i.e. zn[0]) 569 * 570 * cy - c(tn, y) 571 * 572 * c_tmp1, c_tmp2 - temporary N_Vectors (of same length as c) provided 573 * for use by cp_lsetupP. 574 * s_tmp1 - temporary N_Vectors (of same length as y) provided 575 * for use by cp_lsetupP. 576 * 577 * The cp_lsetupP routine should return 0 if successful, a positive 578 * value for a recoverable error, and a negative value for an 579 * unrecoverable error. 580 * ----------------------------------------------------------------- 581 */ 582 583 /* 584 * ----------------------------------------------------------------- 585 * int (*cp_lsolveP)(CPodeMem cp_mem, N_Vector b, N_Vector x, 586 * N_Vector y, N_Vector cy, 587 * N_Vector c_tmp1, N_Vector s_tmp1); 588 * ----------------------------------------------------------------- 589 * cp_lsolveP must solve the linear equation: 590 * 591 * 592 * Explicit-form ODE: 593 * P is some approximation to (I - gamma * df/dy) 594 * and the RHS vector b is input. 595 * Implicit-form ODE: 596 * P is some approximation to (dF/dy' + gamma * dF/dy) 597 * and b is the current residual. 598 * Its arguments are as follows: 599 * 600 * cp_mem - problem memory pointer of type CPodeMem. See the big 601 * typedef earlier in this file. 602 * 603 * b - on input the right-hand side of the linear system 604 * to be solved. On output, the solution. 605 * 606 * weight - te current weights in the WRMS norm. 607 * 608 * yC - the solver's current approximation to y(tn) 609 * 610 * ypC - the solver's current approximation to y'(tn) 611 * 612 * NOTE: ypC = NULL for explicit-form ODE 613 * 614 * fctC - f(tn, yC) or F(tn, yC, ypc). 615 * 616 * The solution is to be returned in the vector b. 617 * The cp_lsolve function should return a positive value for a 618 * recoverable error and a negative value for an unrecoverable error. 619 * Success is indicated by a 0 return value. 620 * ----------------------------------------------------------------- 621 */ 622 623 /* 624 * ----------------------------------------------------------------- 625 * void (*cp_lmultP)(CPodeMem cp_mem, N_Vector x, N_Vector Gx); 626 * ----------------------------------------------------------------- 627 * cp_lmultP should compute the Jacobian - vector product, for a 628 * given vector x and return the result in Gx. 629 * ----------------------------------------------------------------- 630 */ 631 632 /* 633 * ----------------------------------------------------------------- 634 * void (*cp_lfreeP)(CPodeMem cp_mem); 635 * ----------------------------------------------------------------- 636 * cp_lfreeP should free up any memory allocated by the linear 637 * solver. This routine is called once a problem has been 638 * completed and the linear solver is no longer needed. 639 * ----------------------------------------------------------------- 640 */ 641 642 /* 643 * ----------------------------------------------------------------- 644 * Communication between CPODES and a CPODES Linear Solver 645 * ----------------------------------------------------------------- 646 * convfail is passed as an input to cp_lsetup. When dealing with 647 * an ODE in explicit form, this can be used to test whether the Jacobian 648 * must be reevakluated or whether it is enough to use the up-to-date gamma 649 * value. 650 * 651 * convfail can be one of: 652 * 653 * CP_NO_FAILURES : Either this is the first cp_setup call for this 654 * step, or the local error test failed on the 655 * previous attempt at this step (but the Newton 656 * iteration converged). 657 * 658 * CP_FAIL_BAD_J : This value is passed to cp_lsetup if 659 * 660 * (a) The previous Newton corrector iteration 661 * did not converge and the linear solver's 662 * setup routine indicated that its Jacobian- 663 * related data is not current 664 * or 665 * (b) During the previous Newton corrector 666 * iteration, the linear solver's solve routine 667 * failed in a recoverable manner and the 668 * linear solver's setup routine indicated that 669 * its Jacobian-related data is not current. 670 * 671 * CP_FAIL_OTHER : During the current internal step try, the 672 * previous Newton iteration failed to converge 673 * even though the linear solver was using current 674 * Jacobian-related data. 675 * ----------------------------------------------------------------- 676 */ 677 678 679 #define CP_NO_FAILURES 0 680 #define CP_FAIL_BAD_J 1 681 #define CP_FAIL_OTHER 2 682 683 684 #ifdef __cplusplus 685 } 686 #endif 687 688 #endif 689