1 /*-----------------------------------------------------------------
2  * Programmer(s): Daniel R. Reynolds @ SMU
3  *                David J. Gardner, Radu Serban and Aaron Collier @ 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 KINSOL's linear solver interface.
16  *-----------------------------------------------------------------*/
17 
18 #ifndef _KINLS_IMPL_H
19 #define _KINLS_IMPL_H
20 
21 #include <kinsol/kinsol_ls.h>
22 #include "kinsol_impl.h"
23 
24 #ifdef __cplusplus  /* wrapper to enable C++ usage */
25 extern "C" {
26 #endif
27 
28 
29 /*------------------------------------------------------------------
30   keys for KINPrintInfo (do not use 1 -> conflict with PRNT_RETVAL)
31   ------------------------------------------------------------------*/
32 #define PRNT_NLI   101
33 #define PRNT_EPS   102
34 
35 
36 /*------------------------------------------------------------------
37   Types : struct KINLsMemRec, struct *KINLsMem
38 
39   The type KINLsMem is a pointer to a KINLsMemRec, which is a
40   structure containing fields that must be accessible by LS module
41   routines.
42   ------------------------------------------------------------------*/
43 typedef struct KINLsMemRec {
44 
45   /* Linear solver type information */
46   booleantype iterative;    /* is the solver iterative?    */
47   booleantype matrixbased;  /* is a matrix structure used? */
48 
49   /* Jacobian construction & storage */
50   booleantype jacDQ;   /* SUNTRUE if using internal DQ Jacobian approx. */
51   KINLsJacFn jac;      /* Jacobian routine to be called                 */
52   void *J_data;        /* J_data is passed to jac                       */
53 
54   /* Linear solver, matrix and vector objects/pointers */
55   SUNLinearSolver LS;  /* generic iterative linear solver object        */
56   SUNMatrix J;         /* problem Jacobian                              */
57 
58   /* Solver tolerance adjustment factor (if needed, see kinLsSolve)     */
59   realtype tol_fac;
60 
61   /* Statistics and associated parameters */
62   long int nje;        /* no. of calls to jac                           */
63   long int nfeDQ;      /* no. of calls to F due to DQ Jacobian or J*v
64                           approximations                                */
65   long int npe;        /* npe = total number of precond calls           */
66   long int nli;        /* nli = total number of linear iterations       */
67   long int nps;        /* nps = total number of psolve calls            */
68   long int ncfl;       /* ncfl = total number of convergence failures   */
69   long int njtimes;    /* njtimes = total number of calls to jtimes     */
70 
71   booleantype new_uu;  /* flag indicating if the iterate has been
72                           updated - the Jacobian must be updated or
73                           reevaluated (meant to be used by a
74                           user-supplied jtimes function                 */
75 
76   int last_flag;       /* last error return flag                        */
77 
78   /* Preconditioner computation
79      (a) user-provided:
80          - pdata == user_data
81          - pfree == NULL (the user dealocates memory)
82      (b) internal preconditioner module
83          - pdata == kin_mem
84          - pfree == set by the prec. module and called in kinLsFree */
85   KINLsPrecSetupFn pset;
86   KINLsPrecSolveFn psolve;
87   int (*pfree)(KINMem kin_mem);
88   void *pdata;
89 
90   /* Jacobian times vector compuation
91      (a) jtimes function provided by the user:
92          - jt_data == user_data
93          - jtimesDQ == SUNFALSE
94      (b) internal jtimes
95          - jt_data == kin_mem
96          - jtimesDQ == SUNTRUE */
97   booleantype jtimesDQ;
98   KINLsJacTimesVecFn jtimes;
99   KINSysFn jt_func;
100   void *jt_data;
101 
102 } *KINLsMem;
103 
104 
105 /*------------------------------------------------------------------
106   Prototypes of internal functions
107   ------------------------------------------------------------------*/
108 
109 /* Interface routines called by system SUNLinearSolvers */
110 int kinLsATimes(void *kinmem, N_Vector v, N_Vector z);
111 int kinLsPSetup(void *kinmem);
112 int kinLsPSolve(void *kinmem, N_Vector r, N_Vector z,
113                 realtype tol, int lr);
114 
115 /* Difference quotient approximation for Jacobian times vector */
116 int kinLsDQJtimes(N_Vector v, N_Vector Jv, N_Vector u,
117                   booleantype *new_u, void *data);
118 
119 /* Difference-quotient Jacobian approximation routines */
120 int kinLsDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
121                void *data, N_Vector tmp1, N_Vector tmp2);
122 
123 int kinLsDenseDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
124                     KINMem kin_mem, N_Vector tmp1, N_Vector tmp2);
125 
126 int kinLsBandDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
127                    KINMem kin_mem, N_Vector tmp1, N_Vector tmp2);
128 
129 /* Generic linit/lsetup/lsolve/lfree interface routines for KINSOL to call */
130 int kinLsInitialize(KINMem kin_mem);
131 int kinLsSetup(KINMem kin_mem);
132 int kinLsSolve(KINMem kin_mem, N_Vector x, N_Vector b,
133                realtype *sJpnorm, realtype *sFdotJp);
134 int kinLsFree(KINMem kin_mem);
135 
136 /* Auxilliary functions */
137 int kinLsInitializeCounters(KINLsMem kinls_mem);
138 int kinLs_AccessLMem(void* kinmem, const char *fname,
139                      KINMem* kin_mem, KINLsMem *kinls_mem);
140 
141 
142 /*------------------------------------------------------------------
143   Error messages
144   ------------------------------------------------------------------*/
145 
146 #define MSG_LS_KINMEM_NULL    "KINSOL memory is NULL."
147 #define MSG_LS_MEM_FAIL       "A memory request failed."
148 #define MSG_LS_BAD_NVECTOR    "A required vector operation is not implemented."
149 #define MSG_LS_LMEM_NULL      "Linear solver memory is NULL."
150 #define MSG_LS_NEG_MAXRS      "maxrs < 0 illegal."
151 #define MSG_LS_BAD_SIZES      "Illegal bandwidth parameter(s). Must have 0 <=  ml, mu <= N-1."
152 
153 #define MSG_LS_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner."
154 #define MSG_LS_PSET_FAILED    "The preconditioner setup routine failed in an unrecoverable manner."
155 #define MSG_LS_PSOLVE_FAILED  "The preconditioner solve routine failed in an unrecoverable manner."
156 #define MSG_LS_JTIMES_FAILED  "The Jacobian x vector routine failed in an unrecoverable manner."
157 #define MSG_LS_MATZERO_FAILED "The SUNMatZero routine failed in an unrecoverable manner."
158 
159 
160 /*------------------------------------------------------------------
161   Info messages
162   ------------------------------------------------------------------*/
163 
164 #define INFO_NLI  "nli_inc = %d"
165 
166 #if defined(SUNDIALS_EXTENDED_PRECISION)
167 
168 #define INFO_EPS  "residual norm = %12.3Lg  eps = %12.3Lg"
169 
170 #elif defined(SUNDIALS_DOUBLE_PRECISION)
171 
172 #define INFO_EPS  "residual norm = %12.3lg  eps = %12.3lg"
173 
174 #else
175 
176 #define INFO_EPS  "residual norm = %12.3g  eps = %12.3g"
177 
178 #endif
179 
180 
181 #ifdef __cplusplus
182 }
183 #endif
184 
185 #endif
186