1 /*----------------------------------------------------------------- 2 * Programmer(s): Daniel R. Reynolds @ SMU 3 * Alan C. Hindmarsh and Radu Serban @ LLNL 4 *----------------------------------------------------------------- 5 * SUNDIALS Copyright Start 6 * Copyright (c) 2002-2020, Lawrence Livermore National Security 7 * and Southern Methodist University. 8 * All rights reserved. 9 * 10 * See the top-level LICENSE and NOTICE files for details. 11 * 12 * SPDX-License-Identifier: BSD-3-Clause 13 * SUNDIALS Copyright End 14 *----------------------------------------------------------------- 15 * Implementation header file for IDA's linear solver interface. 16 *-----------------------------------------------------------------*/ 17 18 #ifndef _IDALS_IMPL_H 19 #define _IDALS_IMPL_H 20 21 #include <ida/ida_ls.h> 22 #include "ida_impl.h" 23 24 #ifdef __cplusplus /* wrapper to enable C++ usage */ 25 extern "C" { 26 #endif 27 28 /*----------------------------------------------------------------- 29 Types : struct IDALsMemRec, struct *IDALsMem 30 31 The type IDALsMem is a pointer to a IDALsMemRec, which is a 32 structure containing fields that must be accessible by LS module 33 routines. 34 -----------------------------------------------------------------*/ 35 typedef struct IDALsMemRec { 36 37 /* Linear solver type information */ 38 booleantype iterative; /* is the solver iterative? */ 39 booleantype matrixbased; /* is a matrix structure used? */ 40 41 /* Jacobian construction & storage */ 42 booleantype jacDQ; /* SUNTRUE if using internal DQ Jacobian approx. */ 43 IDALsJacFn jac; /* Jacobian routine to be called */ 44 void *J_data; /* J_data is passed to jac */ 45 46 /* Linear solver, matrix and vector objects/pointers */ 47 SUNLinearSolver LS; /* generic linear solver object */ 48 SUNMatrix J; /* J = dF/dy + cj*dF/dy' */ 49 N_Vector ytemp; /* temp vector used by IDAAtimesDQ */ 50 N_Vector yptemp; /* temp vector used by IDAAtimesDQ */ 51 N_Vector x; /* temp vector used by the solve function */ 52 N_Vector ycur; /* current y vector in Newton iteration */ 53 N_Vector ypcur; /* current yp vector in Newton iteration */ 54 N_Vector rcur; /* rcur = F(tn, ycur, ypcur) */ 55 56 /* Matrix-based solver, scale solution to account for change in cj */ 57 booleantype scalesol; 58 59 /* Iterative solver tolerance */ 60 realtype sqrtN; /* sqrt(N) */ 61 realtype eplifac; /* eplifac = linear convergence factor */ 62 63 /* Statistics and associated parameters */ 64 realtype dqincfac; /* dqincfac = optional increment factor in Jv */ 65 long int nje; /* nje = no. of calls to jac */ 66 long int npe; /* npe = total number of precond calls */ 67 long int nli; /* nli = total number of linear iterations */ 68 long int nps; /* nps = total number of psolve calls */ 69 long int ncfl; /* ncfl = total number of convergence failures */ 70 long int nreDQ; /* nreDQ = total number of calls to res */ 71 long int njtsetup; /* njtsetup = total number of calls to jtsetup */ 72 long int njtimes; /* njtimes = total number of calls to jtimes */ 73 long int nst0; /* nst0 = saved nst (for performance monitor) */ 74 long int nni0; /* nni0 = saved nni (for performance monitor) */ 75 long int ncfn0; /* ncfn0 = saved ncfn (for performance monitor) */ 76 long int ncfl0; /* ncfl0 = saved ncfl (for performance monitor) */ 77 long int nwarn; /* nwarn = no. of warnings (for perf. monitor) */ 78 79 int last_flag; /* last error return flag */ 80 81 /* Preconditioner computation 82 (a) user-provided: 83 - pdata == user_data 84 - pfree == NULL (the user dealocates memory) 85 (b) internal preconditioner module 86 - pdata == ida_mem 87 - pfree == set by the prec. module and called in idaLsFree */ 88 IDALsPrecSetupFn pset; 89 IDALsPrecSolveFn psolve; 90 int (*pfree)(IDAMem IDA_mem); 91 void *pdata; 92 93 /* Jacobian times vector compuation 94 (a) jtimes function provided by the user: 95 - jt_data == user_data 96 - jtimesDQ == SUNFALSE 97 (b) internal jtimes 98 - jt_data == ida_mem 99 - jtimesDQ == SUNTRUE */ 100 booleantype jtimesDQ; 101 IDALsJacTimesSetupFn jtsetup; 102 IDALsJacTimesVecFn jtimes; 103 void *jt_data; 104 105 } *IDALsMem; 106 107 108 /*----------------------------------------------------------------- 109 Prototypes of internal functions 110 -----------------------------------------------------------------*/ 111 112 /* Interface routines called by system SUNLinearSolver */ 113 int idaLsATimes(void *ida_mem, N_Vector v, N_Vector z); 114 int idaLsPSetup(void *ida_mem); 115 int idaLsPSolve(void *ida_mem, N_Vector r, N_Vector z, 116 realtype tol, int lr); 117 118 /* Difference quotient approximation for Jac times vector */ 119 int idaLsDQJtimes(realtype tt, N_Vector yy, N_Vector yp, 120 N_Vector rr, N_Vector v, N_Vector Jv, 121 realtype c_j, void *data, 122 N_Vector work1, N_Vector work2); 123 124 /* Difference-quotient Jacobian approximation routines */ 125 int idaLsDQJac(realtype tt, realtype c_j, N_Vector yy, N_Vector yp, 126 N_Vector rr, SUNMatrix Jac, void *data, 127 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 128 int idaLsDenseDQJac(realtype tt, realtype c_j, N_Vector yy, 129 N_Vector yp, N_Vector rr, SUNMatrix Jac, 130 IDAMem IDA_mem, N_Vector tmp1); 131 int idaLsBandDQJac(realtype tt, realtype c_j, N_Vector yy, 132 N_Vector yp, N_Vector rr, SUNMatrix Jac, 133 IDAMem IDA_mem, N_Vector tmp1, 134 N_Vector tmp2, N_Vector tmp3); 135 136 /* Generic linit/lsetup/lsolve/lperf/lfree interface routines for IDA to call */ 137 int idaLsInitialize(IDAMem IDA_mem); 138 int idaLsSetup(IDAMem IDA_mem, N_Vector y, N_Vector yp, N_Vector r, 139 N_Vector vt1, N_Vector vt2, N_Vector vt3); 140 int idaLsSolve(IDAMem IDA_mem, N_Vector b, N_Vector weight, 141 N_Vector ycur, N_Vector ypcur, N_Vector rescur); 142 int idaLsPerf(IDAMem IDA_mem, int perftask); 143 int idaLsFree(IDAMem IDA_mem); 144 145 146 /* Auxilliary functions */ 147 int idaLsInitializeCounters(IDALsMem idals_mem); 148 int idaLs_AccessLMem(void* ida_mem, const char* fname, 149 IDAMem* IDA_mem, IDALsMem* idals_mem); 150 151 152 /*--------------------------------------------------------------- 153 Error and Warning Messages 154 ---------------------------------------------------------------*/ 155 156 #if defined(SUNDIALS_EXTENDED_PRECISION) 157 #define MSG_LS_TIME "at t = %Lg, " 158 #define MSG_LS_FRMT "%Le." 159 #elif defined(SUNDIALS_DOUBLE_PRECISION) 160 #define MSG_LS_TIME "at t = %lg, " 161 #define MSG_LS_FRMT "%le." 162 #else 163 #define MSG_LS_TIME "at t = %g, " 164 #define MSG_LS_FRMT "%e." 165 #endif 166 167 /* Error Messages */ 168 #define MSG_LS_IDAMEM_NULL "Integrator memory is NULL." 169 #define MSG_LS_MEM_FAIL "A memory request failed." 170 #define MSG_LS_BAD_NVECTOR "A required vector operation is not implemented." 171 #define MSG_LS_BAD_SIZES "Illegal bandwidth parameter(s). Must have 0 <= ml, mu <= N-1." 172 #define MSG_LS_BAD_LSTYPE "Incompatible linear solver type." 173 #define MSG_LS_LMEM_NULL "Linear solver memory is NULL." 174 #define MSG_LS_BAD_GSTYPE "gstype has an illegal value." 175 #define MSG_LS_NEG_MAXRS "maxrs < 0 illegal." 176 #define MSG_LS_NEG_EPLIFAC "eplifac < 0.0 illegal." 177 #define MSG_LS_NEG_DQINCFAC "dqincfac < 0.0 illegal." 178 #define MSG_LS_PSET_FAILED "The preconditioner setup routine failed in an unrecoverable manner." 179 #define MSG_LS_PSOLVE_FAILED "The preconditioner solve routine failed in an unrecoverable manner." 180 #define MSG_LS_JTSETUP_FAILED "The Jacobian x vector setup routine failed in an unrecoverable manner." 181 #define MSG_LS_JTIMES_FAILED "The Jacobian x vector routine failed in an unrecoverable manner." 182 #define MSG_LS_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner." 183 #define MSG_LS_MATZERO_FAILED "The SUNMatZero routine failed in an unrecoverable manner." 184 185 /* Warning Messages */ 186 #define MSG_LS_WARN "Warning: " MSG_LS_TIME "poor iterative algorithm performance. " 187 #define MSG_LS_CFN_WARN MSG_LS_WARN "Nonlinear convergence failure rate is " MSG_LS_FRMT 188 #define MSG_LS_CFL_WARN MSG_LS_WARN "Linear convergence failure rate is " MSG_LS_FRMT 189 190 191 #ifdef __cplusplus 192 } 193 #endif 194 195 #endif 196