source: branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/sundials/include/sundials/sundials_spgmr.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: 12.2 KB
Line 
1/*
2 * -----------------------------------------------------------------
3 * $Revision$
4 * $Date$
5 * -----------------------------------------------------------------
6 * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and
7 *                Radu Serban @ LLNL
8 * -----------------------------------------------------------------
9 * LLNS Copyright Start
10 * Copyright (c) 2014, Lawrence Livermore National Security
11 * This work was performed under the auspices of the U.S. Department
12 * of Energy by Lawrence Livermore National Laboratory in part under
13 * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
14 * Produced at the Lawrence Livermore National Laboratory.
15 * All rights reserved.
16 * For details, see the LICENSE file.
17 * LLNS Copyright End
18 * -----------------------------------------------------------------
19 * This is the header file for the implementation of SPGMR Krylov
20 * iterative linear solver.  The SPGMR algorithm is based on the
21 * Scaled Preconditioned GMRES (Generalized Minimal Residual)
22 * method.
23 *
24 * The SPGMR algorithm solves a linear system A x = b.
25 * Preconditioning is allowed on the left, right, or both.
26 * Scaling is allowed on both sides, and restarts are also allowed.
27 * We denote the preconditioner and scaling matrices as follows:
28 *   P1 = left preconditioner
29 *   P2 = right preconditioner
30 *   S1 = diagonal matrix of scale factors for P1-inverse b
31 *   S2 = diagonal matrix of scale factors for P2 x
32 * The matrices A, P1, and P2 are not required explicitly; only
33 * routines that provide A, P1-inverse, and P2-inverse as
34 * operators are required.
35 *
36 * In this notation, SPGMR applies the underlying GMRES method to
37 * the equivalent transformed system
38 *   Abar xbar = bbar , where
39 *   Abar = S1 (P1-inverse) A (P2-inverse) (S2-inverse) ,
40 *   bbar = S1 (P1-inverse) b , and   xbar = S2 P2 x .
41 *
42 * The scaling matrices must be chosen so that vectors S1
43 * P1-inverse b and S2 P2 x have dimensionless components.
44 * If preconditioning is done on the left only (P2 = I), by a
45 * matrix P, then S2 must be a scaling for x, while S1 is a
46 * scaling for P-inverse b, and so may also be taken as a scaling
47 * for x.  Similarly, if preconditioning is done on the right only
48 * (P1 = I, P2 = P), then S1 must be a scaling for b, while S2 is
49 * a scaling for P x, and may also be taken as a scaling for b.
50 *
51 * The stopping test for the SPGMR iterations is on the L2 norm of
52 * the scaled preconditioned residual:
53 *      || bbar - Abar xbar ||_2  <  delta
54 * with an input test constant delta.
55 *
56 * The usage of this SPGMR solver involves supplying two routines
57 * and making three calls.  The user-supplied routines are
58 *    atimes (A_data, x, y) to compute y = A x, given x,
59 * and
60 *    psolve (P_data, y, x, lr)
61 *                to solve P1 x = y or P2 x = y for x, given y.
62 * The three user calls are:
63 *    mem  = SpgmrMalloc(lmax, vec_tmpl);
64 *           to initialize memory,
65 *    flag = SpgmrSolve(mem,A_data,x,b,...,
66 *                      P_data,s1,s2,atimes,psolve,...);
67 *           to solve the system, and
68 *    SpgmrFree(mem);
69 *           to free the memory created by SpgmrMalloc.
70 * Complete details for specifying atimes and psolve and for the
71 * usage calls are given below and in sundials_iterative.h.
72 * -----------------------------------------------------------------
73 */
74
75#ifndef _SPGMR_H
76#define _SPGMR_H
77
78#include <sundials/sundials_iterative.h>
79
80#ifdef __cplusplus  /* wrapper to enable C++ usage */
81extern "C" {
82#endif
83
84/*
85 * -----------------------------------------------------------------
86 * Types: SpgmrMemRec, SpgmrMem
87 * -----------------------------------------------------------------
88 * SpgmrMem is a pointer to an SpgmrMemRec which contains
89 * the memory needed by SpgmrSolve. The SpgmrMalloc routine
90 * returns a pointer of type SpgmrMem which should then be passed
91 * in subsequent calls to SpgmrSolve. The SpgmrFree routine frees
92 * the memory allocated by SpgmrMalloc.
93 *
94 * l_max is the maximum Krylov dimension that SpgmrSolve will be
95 * permitted to use.
96 *
97 * V is the array of Krylov basis vectors v_1, ..., v_(l_max+1),
98 * stored in V[0], ..., V[l_max], where l_max is the second
99 * parameter to SpgmrMalloc. Each v_i is a vector of type
100 * N_Vector.
101 *
102 * Hes is the (l_max+1) x l_max Hessenberg matrix. It is stored
103 * row-wise so that the (i,j)th element is given by Hes[i][j].
104 *
105 * givens is a length 2*l_max array which represents the
106 * Givens rotation matrices that arise in the algorithm. The
107 * Givens rotation matrices F_0, F_1, ..., F_j, where F_i is
108 *
109 *             1
110 *               1
111 *                 c_i  -s_i      <--- row i
112 *                 s_i   c_i
113 *                           1
114 *                             1
115 *
116 * are represented in the givens vector as
117 * givens[0]=c_0, givens[1]=s_0, givens[2]=c_1, givens[3]=s_1,
118 * ..., givens[2j]=c_j, givens[2j+1]=s_j.
119 *
120 * xcor is a vector (type N_Vector) which holds the scaled,
121 * preconditioned correction to the initial guess.
122 *
123 * yg is a length (l_max+1) array of realtype used to hold "short"
124 * vectors (e.g. y and g).
125 *
126 * vtemp is a vector (type N_Vector) used as temporary vector
127 * storage during calculations.
128 * -----------------------------------------------------------------
129 */
130 
131typedef struct _SpgmrMemRec {
132
133  int l_max;
134
135  N_Vector *V;
136  realtype **Hes;
137  realtype *givens;
138  N_Vector xcor;
139  realtype *yg;
140  N_Vector vtemp;
141
142} SpgmrMemRec, *SpgmrMem;
143
144/*
145 * -----------------------------------------------------------------
146 * Function : SpgmrMalloc
147 * -----------------------------------------------------------------
148 * SpgmrMalloc allocates the memory used by SpgmrSolve. It
149 * returns a pointer of type SpgmrMem which the user of the
150 * SPGMR package should pass to SpgmrSolve. The parameter l_max
151 * is the maximum Krylov dimension that SpgmrSolve will be
152 * permitted to use. The parameter vec_tmpl is a pointer to an
153 * N_Vector used as a template to create new vectors by duplication.
154 * This routine returns NULL if there is a memory request failure.
155 * -----------------------------------------------------------------
156 */
157
158SUNDIALS_EXPORT SpgmrMem SpgmrMalloc(int l_max, N_Vector vec_tmpl);
159
160/*
161 * -----------------------------------------------------------------
162 * Function : SpgmrSolve
163 * -----------------------------------------------------------------
164 * SpgmrSolve solves the linear system Ax = b using the SPGMR
165 * method. The return values are given by the symbolic constants
166 * below. The first SpgmrSolve parameter is a pointer to memory
167 * allocated by a prior call to SpgmrMalloc.
168 *
169 * mem is the pointer returned by SpgmrMalloc to the structure
170 * containing the memory needed by SpgmrSolve.
171 *
172 * A_data is a pointer to information about the coefficient
173 * matrix A. This pointer is passed to the user-supplied function
174 * atimes.
175 *
176 * x is the initial guess x_0 upon entry and the solution
177 * N_Vector upon exit with return value SPGMR_SUCCESS or
178 * SPGMR_RES_REDUCED. For all other return values, the output x
179 * is undefined.
180 *
181 * b is the right hand side N_Vector. It is undisturbed by this
182 * function.
183 *
184 * pretype is the type of preconditioning to be used. Its
185 * legal values are enumerated in sundials_iterative.h. These
186 * values are PREC_NONE=0, PREC_LEFT=1, PREC_RIGHT=2, and
187 * PREC_BOTH=3.
188 *
189 * gstype is the type of Gram-Schmidt orthogonalization to be
190 * used. Its legal values are enumerated in sundials_iterative.h.
191 * These values are MODIFIED_GS=0 and CLASSICAL_GS=1.
192 *
193 * delta is the tolerance on the L2 norm of the scaled,
194 * preconditioned residual. On return with value SPGMR_SUCCESS,
195 * this residual satisfies || s1 P1_inv (b - Ax) ||_2 <= delta.
196 *
197 * max_restarts is the maximum number of times the algorithm is
198 * allowed to restart.
199 *
200 * P_data is a pointer to preconditioner information. This
201 * pointer is passed to the user-supplied function psolve.
202 *
203 * s1 is an N_Vector of positive scale factors for P1-inv b, where
204 * P1 is the left preconditioner. (Not tested for positivity.)
205 * Pass NULL if no scaling on P1-inv b is required.
206 *
207 * s2 is an N_Vector of positive scale factors for P2 x, where
208 * P2 is the right preconditioner. (Not tested for positivity.)
209 * Pass NULL if no scaling on P2 x is required.
210 *
211 * atimes is the user-supplied function which performs the
212 * operation of multiplying A by a given vector. Its description
213 * is given in sundials_iterative.h.
214 *
215 * psolve is the user-supplied function which solves a
216 * preconditioner system Pz = r, where P is P1 or P2. Its full
217 * description is given in sundials_iterative.h. The psolve function
218 * will not be called if pretype is NONE; in that case, the user
219 * should pass NULL for psolve.
220 *
221 * res_norm is a pointer to the L2 norm of the scaled,
222 * preconditioned residual. On return with value SPGMR_SUCCESS or
223 * SPGMR_RES_REDUCED, (*res_norm) contains the value
224 * || s1 P1_inv (b - Ax) ||_2 for the computed solution x.
225 * For all other return values, (*res_norm) is undefined. The
226 * caller is responsible for allocating the memory (*res_norm)
227 * to be filled in by SpgmrSolve.
228 *
229 * nli is a pointer to the number of linear iterations done in
230 * the execution of SpgmrSolve. The caller is responsible for
231 * allocating the memory (*nli) to be filled in by SpgmrSolve.
232 *
233 * nps is a pointer to the number of calls made to psolve during
234 * the execution of SpgmrSolve. The caller is responsible for
235 * allocating the memory (*nps) to be filled in by SpgmrSolve.
236 *
237 * Note: Repeated calls can be made to SpgmrSolve with varying
238 * input arguments. If, however, the problem size N or the
239 * maximum Krylov dimension l_max changes, then a call to
240 * SpgmrMalloc must be made to obtain new memory for SpgmrSolve
241 * to use.
242 * -----------------------------------------------------------------
243 */                                                               
244     
245SUNDIALS_EXPORT int SpgmrSolve(SpgmrMem mem, void *A_data, N_Vector x, N_Vector b,
246             int pretype, int gstype, realtype delta,
247             int max_restarts, void *P_data, N_Vector s1,
248             N_Vector s2, ATimesFn atimes, PSolveFn psolve,
249             realtype *res_norm, int *nli, int *nps);
250
251
252/* Return values for SpgmrSolve */
253
254#define SPGMR_SUCCESS            0  /* Converged                     */
255#define SPGMR_RES_REDUCED        1  /* Did not converge, but reduced
256                                       norm of residual              */
257#define SPGMR_CONV_FAIL          2  /* Failed to converge            */
258#define SPGMR_QRFACT_FAIL        3  /* QRfact found singular matrix  */
259#define SPGMR_PSOLVE_FAIL_REC    4  /* psolve failed recoverably     */
260#define SPGMR_ATIMES_FAIL_REC    5  /* atimes failed recoverably     */
261#define SPGMR_PSET_FAIL_REC      6  /* pset faild recoverably        */
262
263#define SPGMR_MEM_NULL          -1  /* mem argument is NULL          */
264#define SPGMR_ATIMES_FAIL_UNREC -2  /* atimes returned failure flag  */
265#define SPGMR_PSOLVE_FAIL_UNREC -3  /* psolve failed unrecoverably   */
266#define SPGMR_GS_FAIL           -4  /* Gram-Schmidt routine faiuled  */       
267#define SPGMR_QRSOL_FAIL        -5  /* QRsol found singular R        */
268#define SPGMR_PSET_FAIL_UNREC   -6  /* pset failed unrecoverably     */
269
270/*
271 * -----------------------------------------------------------------
272 * Function : SpgmrFree
273 * -----------------------------------------------------------------
274 * SpgmrMalloc frees the memory allocated by SpgmrMalloc. It is
275 * illegal to use the pointer mem after a call to SpgmrFree.
276 * -----------------------------------------------------------------
277 */                                                               
278
279SUNDIALS_EXPORT void SpgmrFree(SpgmrMem mem);
280
281/*
282 * -----------------------------------------------------------------
283 * Macro: SPGMR_VTEMP
284 * -----------------------------------------------------------------
285 * This macro provides access to the work vector vtemp in the
286 * memory block of the SPGMR module.  The argument mem is the
287 * memory pointer returned by SpgmrMalloc, of type SpgmrMem,
288 * and the macro value is of type N_Vector.
289 * On a return from SpgmrSolve with *nli = 0, this vector
290 * contains the scaled preconditioned initial residual,
291 * s1 * P1_inverse * (b - A x_0).
292 * -----------------------------------------------------------------
293 */
294
295#define SPGMR_VTEMP(mem) (mem->vtemp)
296
297#ifdef __cplusplus
298}
299#endif
300
301#endif
Note: See TracBrowser for help on using the repository browser.