forked from Ermentrout/xppaut
-
Notifications
You must be signed in to change notification settings - Fork 1
/
cvspgmr.h
executable file
·322 lines (299 loc) · 19.4 KB
/
cvspgmr.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
/******************************************************************
* *
* File : cvspgmr.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Last Modified : 1 September 1994 *
*----------------------------------------------------------------*
* This is the header file for the CVODE scaled, preconditioned *
* GMRES linear solver, CVSPGMR. *
* *
* Note: The type integer must be large enough to store the value *
* of the linear system size N. *
* *
******************************************************************/
#ifndef _cvspgmr_h
#define _cvspgmr_h
#include <stdio.h>
#include "cvode.h"
#include "spgmr.h"
#include "llnltyps.h"
#include "vector.h"
/******************************************************************
* *
* CVSPGMR solver statistics indices *
*----------------------------------------------------------------*
* The following enumeration gives a symbolic name to each *
* CVSPGMR statistic. The symbolic names are used as indices into *
* the iopt and ropt arrays passed to CVodeMalloc. *
* The CVSPGMR statistics are: *
* *
* iopt[SPGMR_NPE] : number of preconditioner evaluations, *
* i.e. of calls made to user's precond *
* function with jok == FALSE. *
* *
* iopt[SPGMR_NLI] : number of linear iterations. *
* *
* iopt[SPGMR_NPS] : number of calls made to user's psolve *
* function. *
* *
* iopt[SPGMR_NCFL] : number of linear convergence failures. *
* *
* iopt[SPGMR_LRW] : size (in real words) of real workspace *
* vectors and small matrices used by this *
* solver. *
* *
* iopt[SPGMR_LIW] : size (in integer words) of integer *
* workspace vectors used by this solver. *
* *
******************************************************************/
enum { SPGMR_NPE = CVODE_IOPT_SIZE,
SPGMR_NLI, SPGMR_NPS, SPGMR_NCFL, SPGMR_LRW, SPGMR_LIW };
/******************************************************************
* *
* CVSPGMR solver constants *
*----------------------------------------------------------------*
* CVSPGMR_MAXL : default value for the maximum Krylov *
* dimension is MIN(N, CVSPGMR_MAXL) *
* *
* CVSPGMR_MSBPRE : maximum number of steps between *
* preconditioner evaluations *
* *
* CVSPGMR_DGMAX : maximum change in gamma between *
* preconditioner evaluations *
* *
* CVSPGMR_DELT : default value for factor by which the *
* tolerance on the nonlinear iteration is *
* multiplied to get a tolerance on the linear *
* iteration *
* *
******************************************************************/
#define CVSPGMR_MAXL 5
#define CVSPGMR_MSBPRE 50
#define CVSPGMR_DGMAX RCONST(0.2)
#define CVSPGMR_DELT RCONST(0.05)
/******************************************************************
* *
* Type : CVSpgmrPrecondFn *
*----------------------------------------------------------------*
* The user-supplied preconditioner setup function Precond and *
* the user-supplied preconditioner solve function PSolve *
* together must define left and right preconditoner matrices *
* P1 and P2 (either of which may be trivial), such that the *
* product P1*P2 is an approximation to the Newton matrix *
* M = I - gamma*J. Here J is the system Jacobian J = df/dy, *
* and gamma is a scalar proportional to the integration step *
* size h. The solution of systems P z = r, with P = P1 or P2, *
* is to be carried out by the PSolve function, and Precond is *
* to do any necessary setup operations. *
* *
* The user-supplied preconditioner setup function Precond *
* is to evaluate and preprocess any Jacobian-related data *
* needed by the preconditioner solve function PSolve. *
* This might include forming a crude approximate Jacobian, *
* and performing an LU factorization on the resulting *
* approximation to M. This function will not be called in *
* advance of every call to PSolve, but instead will be called *
* only as often as necessary to achieve convergence within the *
* Newton iteration in CVODE. If the PSolve function needs no *
* preparation, the Precond function can be NULL. *
* *
* For greater efficiency, the Precond function may save *
* Jacobian-related data and reuse it, rather than generating it *
* from scratch. In this case, it should use the input flag jok *
* to decide whether to recompute the data, and set the output *
* flag *jcurPtr accordingly. *
* *
* The error weight vector ewt, step size h, and unit roundoff *
* uround are provided to the Precond function for possible use *
* in approximating Jacobian data, e.g. by difference quotients. *
* *
* A function Precond must have the prototype given below. *
* Its parameters are as follows: *
* *
* N is the length of all vector arguments. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector, *
* namely the predicted value of y(t). *
* *
* fy is the vector f(t,y). *
* *
* jok is an input flag indicating whether Jacobian-related *
* data needs to be recomputed, as follows: *
* jok == FALSE means recompute Jacobian-related data *
* from scratch. *
* jok == TRUE means that Jacobian data, if saved from *
* the previous Precond call, can be reused *
* (with the current value of gamma). *
* A Precond call with jok == TRUE can only occur after *
* a call with jok == FALSE. *
* *
* jcurPtr is a pointer to an output integer flag which is *
* to be set by Precond as follows: *
* Set *jcurPtr = TRUE if Jacobian data was recomputed. *
* Set *jcurPtr = FALSE if Jacobian data was not *
* recomputed, but saved data was reused. *
* *
* gamma is the scalar appearing in the Newton matrix. *
* *
* ewt is the error weight vector. *
* *
* h is a tentative step size in t. *
* *
* uround is the machine unit roundoff. *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. *
* The Precond routine should update this counter by *
* adding on the number of f calls made in order to *
* approximate the Jacobian, if any. For example, if *
* the routine calls f a total of W times, then the *
* update is *nfePtr += W. *
* *
* P_data is a pointer to user data - the same as the P_data *
* parameter passed to CVSpgmr. *
* *
* vtemp1, vtemp2, and vtemp3 are pointers to memory allocated *
* for vectors of length N which can be used by *
* CVSpgmrPrecondFn as temporary storage or work space. *
* *
* *
* Returned value: *
* The value to be returned by the Precond function is a flag *
* indicating whether it was successful. This value should be *
* 0 if successful, *
* > 0 for a recoverable error (step will be retried), *
* < 0 for an unrecoverable error (integration is halted). *
* *
******************************************************************/
typedef int (*CVSpgmrPrecondFn)(integer N, real t, N_Vector y, N_Vector fy,
bool jok, bool *jcurPtr, real gamma,
N_Vector ewt, real h, real uround,
int *nfePtr, void *P_data,
N_Vector vtemp1, N_Vector vtemp2,
N_Vector vtemp3);
/******************************************************************
* *
* Type : CVSpgmrPSolveFn *
*----------------------------------------------------------------*
* The user-supplied preconditioner solve function PSolve *
* is to solve a linear system P z = r in which the matrix P is *
* one of the preconditioner matrices P1 or P2, depending on the *
* type of preconditioning chosen. *
* *
* A function PSolve must have the prototype given below. *
* Its parameters are as follows: *
* *
* N is the length of all vector arguments. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector. *
* *
* fy is the vector f(t,y). *
* *
* vtemp is a pointer to memory allocated for a vector of *
* length N which can be used by PSolve for work space. *
* *
* gamma is the scalar appearing in the Newton matrix. *
* *
* ewt is the error weight vector (input). See delta below. *
* *
* delta is an input tolerance for use by PSolve if it uses *
* an iterative method in its solution. In that case, *
* the residual vector Res = r - P z of the system *
* should be made less than delta in weighted L2 norm, *
* i.e., sqrt [ Sum (Res[i]*ewt[i])^2 ] < delta . *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. The *
* PSolve routine should update this counter by adding *
* on the number of f calls made in order to carry out *
* the solution, if any. For example, if the routine *
* calls f a total of W times, then the update is *
* *nfePtr += W. *
* *
* r is the right-hand side vector of the linear system. *
* *
* lr is an input flag indicating whether PSolve is to use *
* the left preconditioner P1 or right preconditioner *
* P2: lr = 1 means use P1, and lr = 2 means use P2. *
* *
* P_data is a pointer to user data - the same as the P_data *
* parameter passed to CVSpgmr. *
* *
* z is the output vector computed by PSolve. *
* *
* Returned value: *
* The value to be returned by the PSolve function is a flag *
* indicating whether it was successful. This value should be *
* 0 if successful, *
* positive for a recoverable error (step will be retried), *
* negative for an unrecoverable error (integration is halted). *
* *
******************************************************************/
typedef int (*CVSpgmrPSolveFn)(integer N, real t, N_Vector y, N_Vector fy,
N_Vector vtemp, real gamma, N_Vector ewt,
real delta, int *nfePtr, N_Vector r,
int lr, void *P_data, N_Vector z);
/******************************************************************
* *
* Function : CVSpgmr *
*----------------------------------------------------------------*
* A call to the CVSpgmr function links the main CVODE integrator *
* with the CVSPGMR linear solver. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
* pretype is the type of user preconditioning to be done. *
* This must be one of the four enumeration constants *
* NONE, LEFT, RIGHT, or BOTH defined in iterativ.h. *
* These correspond to no preconditioning, *
* left preconditioning only, right preconditioning *
* only, and both left and right preconditioning, *
* respectively. *
* *
* gstype is the type of Gram-Schmidt orthogonalization to be *
* used. This must be one of the two enumeration *
* constants MODIFIED_GS or CLASSICAL_GS defined in *
* iterativ.h. These correspond to using modified *
* Gram-Schmidt and classical Gram-Schmidt, *
* respectively. *
* *
* maxl is the maximum Krylov dimension. This is an *
* optional input to the CVSPGMR solver. Pass 0 to *
* use the default value MIN(N, CVSPGMR_MAXL=5). *
* *
* delt is the factor by which the tolerance on the *
* nonlinear iteration is multiplied to get a *
* tolerance on the linear iteration. This is an *
* optional input to the CVSPGMR solver. Pass 0 to *
* use the default value CVSPGMR_DELT = 0.05. *
* *
* precond is the user's preconditioner routine. It is used to *
* evaluate and preprocess any Jacobian-related data *
* needed by the psolve routine. See the *
* documentation for the type CVSpgmrPrecondFn for *
* full details. Pass NULL if no such setup of *
* Jacobian data is required. A precond routine is *
* NOT required for any of the four possible values *
* of pretype. *
* *
* psolve is the user's preconditioner solve routine. It is *
* used to solve Pz=r, where P is a preconditioner *
* matrix. See the documentation for the type *
* CVSpgmrPSolveFn for full details. The only case *
* in which psolve is allowed to be NULL is when *
* pretype is NONE. A valid psolve function must be *
* supplied when any preconditioning is to be done. *
* *
* P_data is a pointer to user preconditioner data. This *
* pointer is passed to precond and psolve every time *
* these routines are called. *
* *
******************************************************************/
void CVSpgmr(void *cvode_mem, int pretype, int gstype, int maxl, real delt,
CVSpgmrPrecondFn precond, CVSpgmrPSolveFn psolve, void *P_data);
#endif