source: branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/sundials/include/sundials/sundials_spfgmr.h @ 16222

Last change on this file since 16222 was 16222, checked in by gkronber, 4 years ago

#2925:

  • added comments about parameter identification for differential equation models
  • added source code of cvodes library (part of sundials) which provides functionality to calculate gradients for the parameters of partial differential equation models efficiently using the 'adjoint state method'.
  • added compiled version of cvodes
File size: 11.6 KB
Line 
1/* -----------------------------------------------------------------
2 Programmer(s): Daniel R. Reynolds and Hilari C. Tiedeman @ SMU
3 -------------------------------------------------------------------
4 Copyright (c) 2011, Southern Methodist University.
5 All rights reserved.
6 For details, see the LICENSE file.
7 -------------------------------------------------------------------
8 This is the header file for the implementation of SPFGMR Krylov
9 iterative linear solver.  The SPFGMR algorithm is based on the
10 Scaled Preconditioned Flexible GMRES (Generalized Minimal Residual)
11 method [Y. Saad, SIAM J. Sci. Comput., 1993].
12
13 The SPFGMR algorithm solves a linear system A x = b.
14 Preconditioning is only allowed on the right.
15 Scaling is allowed on the right, and restarts are also allowed.
16 We denote the preconditioner and scaling matrices as follows:
17   P = right preconditioner
18   S1 = diagonal matrix of scale factors for P-inverse b
19   S2 = diagonal matrix of scale factors for x
20 The matrices A and P are not required explicitly; only
21 routines that provide A and P-inverse as operators are required.
22
23 In this notation, SPFGMR applies the underlying FGMRES method to
24 the equivalent transformed system
25   Abar xbar = bbar , where
26   Abar = S1 A (P-inverse) (S2-inverse),
27   bbar = S1 b , and   xbar = S2 P x .
28
29 The scaling matrix must be chosen so that the vectors S1 b and
30 S2 P x have dimensionless components.  If preconditioning is not
31 performed (P = I), then S2 must be a scaling for x, while S1 is a
32 scaling for b.  Similarly, if preconditioning is performed, then S1
33 must be a scaling for b, while S2 is a scaling for P x, and may
34 also be taken as a scaling for b.
35
36 The stopping test for the SPFGMR iterations is on the L2 norm of
37 the scaled preconditioned residual:
38      || bbar - Abar xbar ||_2  <  delta
39 with an input test constant delta.
40
41 The usage of this SPFGMR solver involves supplying two routines
42 and making three calls.  The user-supplied routines are
43    atimes (A_data, x, y) to compute y = A x, given x,
44 and
45    psolve (P_data, y, x, lr)
46                to solve P x = y for x, given y.
47 The three user calls are:
48    mem  = SpfgmrMalloc(lmax, vec_tmpl);
49           to initialize memory,
50    flag = SpfgmrSolve(mem,A_data,x,b,...,
51                       P_data,s1,s2,atimes,psolve,...);
52           to solve the system, and
53    SpfgmrFree(mem);
54           to free the memory created by SpfgmrMalloc.
55 Complete details for specifying atimes and psolve and for the
56 usage calls are given in the paragraphs below and in iterative.h.
57 -----------------------------------------------------------------*/
58
59#ifndef _SPFGMR_H
60#define _SPFGMR_H
61
62#include <sundials/sundials_iterative.h>
63
64#ifdef __cplusplus  /* wrapper to enable C++ usage */
65extern "C" {
66#endif
67
68/* -----------------------------------------------------------------
69 Types: SpfgmrMemRec, SpfgmrMem
70 -------------------------------------------------------------------
71 SpfgmrMem is a pointer to an SpfgmrMemRec which contains
72 the memory needed by SpfgmrSolve. The SpfgmrMalloc routine
73 returns a pointer of type SpfgmrMem which should then be passed
74 in subsequent calls to SpfgmrSolve. The SpfgmrFree routine frees
75 the memory allocated by SpfgmrMalloc.
76
77 l_max is the maximum Krylov dimension that SpfgmrSolve will be
78 permitted to use.
79
80 V is the array of Krylov basis vectors v_1, ..., v_(l_max+1),
81 stored in V[0], ..., V[l_max], where l_max is the second
82 parameter to SpfgmrMalloc. Each v_i is a vector of type
83 N_Vector.
84
85 Z is the array of preconditioned basis vectors z_1, ...,
86 z_(l_max+1), stored in Z[0], ..., Z[l_max], where l_max is the
87 second parameter to SpfgmrMalloc. Each z_i is a vector of type
88 N_Vector.
89
90 Hes is the (l_max+1) x l_max Hessenberg matrix. It is stored
91 row-wise so that the (i,j)th element is given by Hes[i][j].
92
93 givens is a length 2*l_max array which represents the
94 Givens rotation matrices that arise in the algorithm. The
95 Givens rotation matrices F_0, F_1, ..., F_j, where F_i is
96
97             1
98               1
99                 c_i  -s_i      <--- row i
100                 s_i   c_i
101                           1
102                             1
103
104 are represented in the givens vector as
105 givens[0]=c_0, givens[1]=s_0, givens[2]=c_1, givens[3]=s_1,
106 ..., givens[2j]=c_j, givens[2j+1]=s_j.
107
108 xcor is a vector (type N_Vector) which holds the scaled,
109 preconditioned correction to the initial guess.
110
111 yg is a length (l_max+1) array of realtype used to hold "short"
112 vectors (e.g. y and g).
113
114 vtemp is a vector (type N_Vector) used as temporary vector
115 storage during calculations.
116 -----------------------------------------------------------------*/
117typedef struct _SpfgmrMemRec {
118  int l_max;
119  N_Vector *V;
120  N_Vector *Z;
121  realtype **Hes;
122  realtype *givens;
123  N_Vector xcor;
124  realtype *yg;
125  N_Vector vtemp;
126} SpfgmrMemRec, *SpfgmrMem;
127
128/*----------------------------------------------------------------
129 Function : SpfgmrMalloc
130 -----------------------------------------------------------------
131 SpfgmrMalloc allocates the memory used by SpfgmrSolve. It
132 returns a pointer of type SpfgmrMem which the user of the
133 SPGMR package should pass to SpfgmrSolve. The parameter l_max
134 is the maximum Krylov dimension that SpfgmrSolve will be
135 permitted to use. The parameter vec_tmpl is a pointer to an
136 N_Vector used as a template to create new vectors by duplication.
137 This routine returns NULL if there is a memory request failure.
138 ---------------------------------------------------------------*/
139
140SUNDIALS_EXPORT SpfgmrMem SpfgmrMalloc(int l_max, N_Vector vec_tmpl);
141
142/*----------------------------------------------------------------
143 Function : SpfgmrSolve
144 -----------------------------------------------------------------
145 SpfgmrSolve solves the linear system Ax = b using the SPFGMR
146 method. The return values are given by the symbolic constants
147 below. The first SpfgmrSolve parameter is a pointer to memory
148 allocated by a prior call to SpfgmrMalloc.
149
150 mem is the pointer returned by SpfgmrMalloc to the structure
151 containing the memory needed by SpfgmrSolve.
152
153 A_data is a pointer to information about the coefficient
154 matrix A. This pointer is passed to the user-supplied function
155 atimes.
156
157 x is the initial guess x_0 upon entry and the solution
158 N_Vector upon exit with return value SPFGMR_SUCCESS or
159 SPFGMR_RES_REDUCED. For all other return values, the output x
160 is undefined.
161
162 b is the right hand side N_Vector. It is undisturbed by this
163 function.
164
165 pretype is the type of preconditioning to be used. Its
166 legal possible values are enumerated in iterative.h. These
167 values are PREC_NONE, PREC_LEFT, PREC_RIGHT and PREC_BOTH;
168 however since this solver can only precondition on the right,
169 then right-preconditioning will be done if any of the values
170 PREC_LEFT, PREC_RIGHT or PREC_BOTH are provided..
171
172 gstype is the type of Gram-Schmidt orthogonalization to be
173 used. Its legal values are enumerated in iterativ.h. These
174 values are MODIFIED_GS=0 and CLASSICAL_GS=1.
175
176 delta is the tolerance on the L2 norm of the scaled,
177 preconditioned residual. On return with value SPFGMR_SUCCESS,
178 this residual satisfies || s1 P1_inv (b - Ax) ||_2 <= delta.
179
180 max_restarts is the maximum number of times the algorithm is
181 allowed to restart.
182
183 maxit is the maximum number of iterations allowed within the
184 solve.  This value must be less than or equal to the "l_max"
185 value previously supplied to SpfgmrMalloc.  If maxit is too
186 large, l_max will be used instead.
187
188 P_data is a pointer to preconditioner information. This
189 pointer is passed to the user-supplied function psolve.
190
191 s1 is an N_Vector of positive scale factors for b. (Not
192 tested for positivity.)  Pass NULL if no scaling on b is
193 required.
194
195 s2 is an N_Vector of positive scale factors for P x, where
196 P is the right preconditioner. (Not tested for positivity.)
197 Pass NULL if no scaling on P x is required.
198
199 atimes is the user-supplied function which performs the
200 operation of multiplying A by a given vector. Its description
201 is given in iterative.h.
202
203 psolve is the user-supplied function which solves a
204 preconditioner system Pz = r, where P is P1 or P2. Its full
205 description is given in iterative.h. The psolve function will
206 not be called if pretype is NONE; in that case, the user
207 should pass NULL for psolve.
208
209 res_norm is a pointer to the L2 norm of the scaled,
210 preconditioned residual. On return with value SPFGMR_SUCCESS or
211 SPFGMR_RES_REDUCED, (*res_norm) contains the value
212 || s1 (b - Ax) ||_2 for the computed solution x.
213 For all other return values, (*res_norm) is undefined. The
214 caller is responsible for allocating the memory (*res_norm)
215 to be filled in by SpfgmrSolve.
216
217 nli is a pointer to the number of linear iterations done in
218 the execution of SpfgmrSolve. The caller is responsible for
219 allocating the memory (*nli) to be filled in by SpfgmrSolve.
220
221 nps is a pointer to the number of calls made to psolve during
222 the execution of SpfgmrSolve. The caller is responsible for
223 allocating the memory (*nps) to be filled in by SpfgmrSolve.
224
225 Note: Repeated calls can be made to SpfgmrSolve with varying
226 input arguments. If, however, the problem size N or the
227 maximum Krylov dimension l_max changes, then a call to
228 SpfgmrMalloc must be made to obtain new memory for SpfgmrSolve
229 to use.
230 ---------------------------------------------------------------*/
231     
232SUNDIALS_EXPORT int SpfgmrSolve(SpfgmrMem mem, void *A_data, N_Vector x,
233        N_Vector b, int pretype, int gstype,
234        realtype delta, int max_restarts,
235        int maxit, void *P_data, N_Vector s1,
236        N_Vector s2, ATimesFn atimes, PSolveFn psolve,
237        realtype *res_norm, int *nli, int *nps);
238
239
240/* Return values for SpfgmrSolve */
241#define SPFGMR_SUCCESS            0  /* Converged                     */
242#define SPFGMR_RES_REDUCED        1  /* Did not converge, but reduced
243                                        norm of residual              */
244#define SPFGMR_CONV_FAIL          2  /* Failed to converge            */
245#define SPFGMR_QRFACT_FAIL        3  /* QRfact found singular matrix  */
246#define SPFGMR_PSOLVE_FAIL_REC    4  /* psolve failed recoverably     */
247#define SPFGMR_ATIMES_FAIL_REC    5  /* atimes failed recoverably     */
248#define SPFGMR_PSET_FAIL_REC      6  /* pset faild recoverably        */
249
250#define SPFGMR_MEM_NULL          -1  /* mem argument is NULL          */
251#define SPFGMR_ATIMES_FAIL_UNREC -2  /* atimes returned failure flag  */
252#define SPFGMR_PSOLVE_FAIL_UNREC -3  /* psolve failed unrecoverably   */
253#define SPFGMR_GS_FAIL           -4  /* Gram-Schmidt routine faiuled  */       
254#define SPFGMR_QRSOL_FAIL        -5  /* QRsol found singular R        */
255#define SPFGMR_PSET_FAIL_UNREC   -6  /* pset failed unrecoverably     */
256
257/*----------------------------------------------------------------
258 Function : SpfgmrFree
259 -----------------------------------------------------------------
260 SpfgmrMalloc frees the memory allocated by SpfgmrMalloc. It is
261 illegal to use the pointer mem after a call to SpfgmrFree.
262 ---------------------------------------------------------------*/
263
264SUNDIALS_EXPORT void SpfgmrFree(SpfgmrMem mem);
265
266/*----------------------------------------------------------------
267 Macro: SPFGMR_VTEMP
268 -----------------------------------------------------------------
269 This macro provides access to the work vector vtemp in the
270 memory block of the SPFGMR module.  The argument mem is the
271 memory pointer returned by SpfgmrMalloc, of type SpfgmrMem,
272 and the macro value is of type N_Vector.
273 On a return from SpfgmrSolve with *nli = 0, this vector
274 contains the scaled preconditioned initial residual,
275 s1 * P1_inverse * (b - A x_0).
276 ---------------------------------------------------------------*/
277
278#define SPFGMR_VTEMP(mem) (mem->vtemp)
279
280#ifdef __cplusplus
281}
282#endif
283
284#endif
Note: See TracBrowser for help on using the repository browser.