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   void *jt_data;
100 
101 } *KINLsMem;
102 
103 
104 /*------------------------------------------------------------------
105   Prototypes of internal functions
106   ------------------------------------------------------------------*/
107 
108 /* Interface routines called by system SUNLinearSolvers */
109 int kinLsATimes(void *kinmem, N_Vector v, N_Vector z);
110 int kinLsPSetup(void *kinmem);
111 int kinLsPSolve(void *kinmem, N_Vector r, N_Vector z,
112                 realtype tol, int lr);
113 
114 /* Difference quotient approximation for Jacobian times vector */
115 int kinLsDQJtimes(N_Vector v, N_Vector Jv, N_Vector u,
116                   booleantype *new_u, void *data);
117 
118 /* Difference-quotient Jacobian approximation routines */
119 int kinLsDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
120                void *data, N_Vector tmp1, N_Vector tmp2);
121 
122 int kinLsDenseDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
123                     KINMem kin_mem, N_Vector tmp1, N_Vector tmp2);
124 
125 int kinLsBandDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
126                    KINMem kin_mem, N_Vector tmp1, N_Vector tmp2);
127 
128 /* Generic linit/lsetup/lsolve/lfree interface routines for KINSOL to call */
129 int kinLsInitialize(KINMem kin_mem);
130 int kinLsSetup(KINMem kin_mem);
131 int kinLsSolve(KINMem kin_mem, N_Vector x, N_Vector b,
132                realtype *sJpnorm, realtype *sFdotJp);
133 int kinLsFree(KINMem kin_mem);
134 
135 /* Auxilliary functions */
136 int kinLsInitializeCounters(KINLsMem kinls_mem);
137 int kinLs_AccessLMem(void* kinmem, const char *fname,
138                      KINMem* kin_mem, KINLsMem *kinls_mem);
139 
140 
141 /*------------------------------------------------------------------
142   Error messages
143   ------------------------------------------------------------------*/
144 
145 #define MSG_LS_KINMEM_NULL    "KINSOL memory is NULL."
146 #define MSG_LS_MEM_FAIL       "A memory request failed."
147 #define MSG_LS_BAD_NVECTOR    "A required vector operation is not implemented."
148 #define MSG_LS_LMEM_NULL      "Linear solver memory is NULL."
149 #define MSG_LS_NEG_MAXRS      "maxrs < 0 illegal."
150 #define MSG_LS_BAD_SIZES      "Illegal bandwidth parameter(s). Must have 0 <=  ml, mu <= N-1."
151 
152 #define MSG_LS_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner."
153 #define MSG_LS_PSET_FAILED    "The preconditioner setup routine failed in an unrecoverable manner."
154 #define MSG_LS_PSOLVE_FAILED  "The preconditioner solve routine failed in an unrecoverable manner."
155 #define MSG_LS_JTIMES_FAILED  "The Jacobian x vector routine failed in an unrecoverable manner."
156 #define MSG_LS_MATZERO_FAILED "The SUNMatZero routine failed in an unrecoverable manner."
157 
158 
159 /*------------------------------------------------------------------
160   Info messages
161   ------------------------------------------------------------------*/
162 
163 #define INFO_NLI  "nli_inc = %d"
164 
165 #if defined(SUNDIALS_EXTENDED_PRECISION)
166 
167 #define INFO_EPS  "residual norm = %12.3Lg  eps = %12.3Lg"
168 
169 #elif defined(SUNDIALS_DOUBLE_PRECISION)
170 
171 #define INFO_EPS  "residual norm = %12.3lg  eps = %12.3lg"
172 
173 #else
174 
175 #define INFO_EPS  "residual norm = %12.3g  eps = %12.3g"
176 
177 #endif
178 
179 
180 #ifdef __cplusplus
181 }
182 #endif
183 
184 #endif
185