1 | /* |
---|
2 | * ----------------------------------------------------------------- |
---|
3 | * $Revision$ |
---|
4 | * $Date$ |
---|
5 | * ----------------------------------------------------------------- |
---|
6 | * Programmer(s): Scott D. Cohen and Alan C. Hindmarsh @ LLNL |
---|
7 | * ----------------------------------------------------------------- |
---|
8 | * LLNS Copyright Start |
---|
9 | * Copyright (c) 2014, Lawrence Livermore National Security |
---|
10 | * This work was performed under the auspices of the U.S. Department |
---|
11 | * of Energy by Lawrence Livermore National Laboratory in part under |
---|
12 | * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344. |
---|
13 | * Produced at the Lawrence Livermore National Laboratory. |
---|
14 | * All rights reserved. |
---|
15 | * For details, see the LICENSE file. |
---|
16 | * LLNS Copyright End |
---|
17 | * ----------------------------------------------------------------- |
---|
18 | * This header file contains declarations intended for use by |
---|
19 | * generic iterative solvers of Ax = b. The enumeration gives |
---|
20 | * symbolic names for the type of preconditioning to be used. |
---|
21 | * The function type declarations give the prototypes for the |
---|
22 | * functions to be called within an iterative linear solver, that |
---|
23 | * are responsible for |
---|
24 | * multiplying A by a given vector v (ATimesFn), |
---|
25 | * setting up a preconditioner P (PSetupFn), and |
---|
26 | * solving the preconditioner equation Pz = r (PSolveFn). |
---|
27 | * ----------------------------------------------------------------- |
---|
28 | */ |
---|
29 | |
---|
30 | #ifndef _ITERATIVE_H |
---|
31 | #define _ITERATIVE_H |
---|
32 | |
---|
33 | #include <sundials/sundials_nvector.h> |
---|
34 | |
---|
35 | #ifdef __cplusplus /* wrapper to enable C++ usage */ |
---|
36 | extern "C" { |
---|
37 | #endif |
---|
38 | |
---|
39 | |
---|
40 | /* |
---|
41 | * ----------------------------------------------------------------- |
---|
42 | * enum : types of preconditioning |
---|
43 | * ----------------------------------------------------------------- |
---|
44 | * PREC_NONE : The iterative linear solver should not use |
---|
45 | * preconditioning. |
---|
46 | * |
---|
47 | * PREC_LEFT : The iterative linear solver uses preconditioning on |
---|
48 | * the left only. |
---|
49 | * |
---|
50 | * PREC_RIGHT : The iterative linear solver uses preconditioning on |
---|
51 | * the right only. |
---|
52 | * |
---|
53 | * PREC_BOTH : The iterative linear solver uses preconditioning on |
---|
54 | * both the left and the right. |
---|
55 | * ----------------------------------------------------------------- |
---|
56 | */ |
---|
57 | |
---|
58 | enum { PREC_NONE, PREC_LEFT, PREC_RIGHT, PREC_BOTH }; |
---|
59 | |
---|
60 | /* |
---|
61 | * ----------------------------------------------------------------- |
---|
62 | * enum : types of Gram-Schmidt routines |
---|
63 | * ----------------------------------------------------------------- |
---|
64 | * MODIFIED_GS : The iterative solver uses the modified |
---|
65 | * Gram-Schmidt routine ModifiedGS listed in this |
---|
66 | * file. |
---|
67 | * |
---|
68 | * CLASSICAL_GS : The iterative solver uses the classical |
---|
69 | * Gram-Schmidt routine ClassicalGS listed in this |
---|
70 | * file. |
---|
71 | * ----------------------------------------------------------------- |
---|
72 | */ |
---|
73 | |
---|
74 | enum { MODIFIED_GS = 1, CLASSICAL_GS = 2 }; |
---|
75 | |
---|
76 | /* |
---|
77 | * ----------------------------------------------------------------- |
---|
78 | * Type: ATimesFn |
---|
79 | * ----------------------------------------------------------------- |
---|
80 | * An ATimesFn multiplies Av and stores the result in z. The |
---|
81 | * caller is responsible for allocating memory for the z vector. |
---|
82 | * The parameter A_data is a pointer to any information about A |
---|
83 | * which the function needs in order to do its job. The vector v |
---|
84 | * is unchanged. An ATimesFn returns 0 if successful and a |
---|
85 | * non-zero value if unsuccessful. |
---|
86 | * ----------------------------------------------------------------- |
---|
87 | */ |
---|
88 | |
---|
89 | typedef int (*ATimesFn)(void *A_data, N_Vector v, N_Vector z); |
---|
90 | |
---|
91 | /* |
---|
92 | * ----------------------------------------------------------------- |
---|
93 | * Type: PSetupFn |
---|
94 | * ----------------------------------------------------------------- |
---|
95 | * A PSetupFn is an integrator-supplied routine that accesses data |
---|
96 | * stored in the integrator memory structure (P_data), and calls |
---|
97 | * the user-supplied, integrator-specific preconditioner setup |
---|
98 | * routine. |
---|
99 | * ----------------------------------------------------------------- |
---|
100 | */ |
---|
101 | |
---|
102 | typedef int (*PSetupFn)(void *P_data); |
---|
103 | |
---|
104 | /* |
---|
105 | * ----------------------------------------------------------------- |
---|
106 | * Type: PSolveFn |
---|
107 | * ----------------------------------------------------------------- |
---|
108 | * A PSolveFn solves the preconditioner equation Pz = r for the |
---|
109 | * vector z. The caller is responsible for allocating memory for |
---|
110 | * the z vector. The parameter P_data is a pointer to any |
---|
111 | * information about P which the function needs in order to do |
---|
112 | * its job. The parameter lr is input, and indicates whether P |
---|
113 | * is to be taken as the left preconditioner or the right |
---|
114 | * preconditioner: lr = 1 for left and lr = 2 for right. |
---|
115 | * If preconditioning is on one side only, lr can be ignored. |
---|
116 | * If the preconditioner is iterative, then it should strive to |
---|
117 | * solve the preconditioner equation so that |
---|
118 | * || Pz - r ||_wrms < tol |
---|
119 | * where the weight vector for the WRMS norm may be accessed from |
---|
120 | * the main integrator memory structure. |
---|
121 | * The vector r should not be modified by the PSolveFn. |
---|
122 | * A PSolveFn returns 0 if successful and a non-zero value if |
---|
123 | * unsuccessful. On a failure, a negative return value indicates |
---|
124 | * an unrecoverable condition, while a positive value indicates |
---|
125 | * a recoverable one, in which the calling routine may reattempt |
---|
126 | * the solution after updating preconditioner data. |
---|
127 | * ----------------------------------------------------------------- |
---|
128 | */ |
---|
129 | |
---|
130 | typedef int (*PSolveFn)(void *P_data, N_Vector r, N_Vector z, |
---|
131 | realtype tol, int lr); |
---|
132 | |
---|
133 | /* |
---|
134 | * ----------------------------------------------------------------- |
---|
135 | * Function: ModifiedGS |
---|
136 | * ----------------------------------------------------------------- |
---|
137 | * ModifiedGS performs a modified Gram-Schmidt orthogonalization |
---|
138 | * of the N_Vector v[k] against the p unit N_Vectors at |
---|
139 | * v[k-1], v[k-2], ..., v[k-p]. |
---|
140 | * |
---|
141 | * v is an array of (k+1) N_Vectors v[i], i=0, 1, ..., k. |
---|
142 | * v[k-1], v[k-2], ..., v[k-p] are assumed to have L2-norm |
---|
143 | * equal to 1. |
---|
144 | * |
---|
145 | * h is the output k by k Hessenberg matrix of inner products. |
---|
146 | * This matrix must be allocated row-wise so that the (i,j)th |
---|
147 | * entry is h[i][j]. The inner products (v[i],v[k]), |
---|
148 | * i=i0, i0+1, ..., k-1, are stored at h[i][k-1]. Here |
---|
149 | * i0=SUNMAX(0,k-p). |
---|
150 | * |
---|
151 | * k is the index of the vector in the v array that needs to be |
---|
152 | * orthogonalized against previous vectors in the v array. |
---|
153 | * |
---|
154 | * p is the number of previous vectors in the v array against |
---|
155 | * which v[k] is to be orthogonalized. |
---|
156 | * |
---|
157 | * new_vk_norm is a pointer to memory allocated by the caller to |
---|
158 | * hold the Euclidean norm of the orthogonalized vector v[k]. |
---|
159 | * |
---|
160 | * If (k-p) < 0, then ModifiedGS uses p=k. The orthogonalized |
---|
161 | * v[k] is NOT normalized and is stored over the old v[k]. Once |
---|
162 | * the orthogonalization has been performed, the Euclidean norm |
---|
163 | * of v[k] is stored in (*new_vk_norm). |
---|
164 | * |
---|
165 | * ModifiedGS returns 0 to indicate success. It cannot fail. |
---|
166 | * ----------------------------------------------------------------- |
---|
167 | */ |
---|
168 | |
---|
169 | SUNDIALS_EXPORT int ModifiedGS(N_Vector *v, realtype **h, int k, int p, |
---|
170 | realtype *new_vk_norm); |
---|
171 | |
---|
172 | /* |
---|
173 | * ----------------------------------------------------------------- |
---|
174 | * Function: ClassicalGS |
---|
175 | * ----------------------------------------------------------------- |
---|
176 | * ClassicalGS performs a classical Gram-Schmidt |
---|
177 | * orthogonalization of the N_Vector v[k] against the p unit |
---|
178 | * N_Vectors at v[k-1], v[k-2], ..., v[k-p]. The parameters v, h, |
---|
179 | * k, p, and new_vk_norm are as described in the documentation |
---|
180 | * for ModifiedGS. |
---|
181 | * |
---|
182 | * temp is an N_Vector which can be used as workspace by the |
---|
183 | * ClassicalGS routine. |
---|
184 | * |
---|
185 | * s is a length k array of realtype which can be used as |
---|
186 | * workspace by the ClassicalGS routine. |
---|
187 | * |
---|
188 | * ClassicalGS returns 0 to indicate success. It cannot fail. |
---|
189 | * ----------------------------------------------------------------- |
---|
190 | */ |
---|
191 | |
---|
192 | SUNDIALS_EXPORT int ClassicalGS(N_Vector *v, realtype **h, int k, int p, |
---|
193 | realtype *new_vk_norm, N_Vector temp, realtype *s); |
---|
194 | |
---|
195 | /* |
---|
196 | * ----------------------------------------------------------------- |
---|
197 | * Function: QRfact |
---|
198 | * ----------------------------------------------------------------- |
---|
199 | * QRfact performs a QR factorization of the Hessenberg matrix H. |
---|
200 | * |
---|
201 | * n is the problem size; the matrix H is (n+1) by n. |
---|
202 | * |
---|
203 | * h is the (n+1) by n Hessenberg matrix H to be factored. It is |
---|
204 | * stored row-wise. |
---|
205 | * |
---|
206 | * q is an array of length 2*n containing the Givens rotations |
---|
207 | * computed by this function. A Givens rotation has the form: |
---|
208 | * | c -s | |
---|
209 | * | s c |. |
---|
210 | * The components of the Givens rotations are stored in q as |
---|
211 | * (c, s, c, s, ..., c, s). |
---|
212 | * |
---|
213 | * job is a control flag. If job==0, then a new QR factorization |
---|
214 | * is performed. If job!=0, then it is assumed that the first |
---|
215 | * n-1 columns of h have already been factored and only the last |
---|
216 | * column needs to be updated. |
---|
217 | * |
---|
218 | * QRfact returns 0 if successful. If a zero is encountered on |
---|
219 | * the diagonal of the triangular factor R, then QRfact returns |
---|
220 | * the equation number of the zero entry, where the equations are |
---|
221 | * numbered from 1, not 0. If QRsol is subsequently called in |
---|
222 | * this situation, it will return an error because it could not |
---|
223 | * divide by the zero diagonal entry. |
---|
224 | * ----------------------------------------------------------------- |
---|
225 | */ |
---|
226 | |
---|
227 | SUNDIALS_EXPORT int QRfact(int n, realtype **h, realtype *q, int job); |
---|
228 | |
---|
229 | /* |
---|
230 | * ----------------------------------------------------------------- |
---|
231 | * Function: QRsol |
---|
232 | * ----------------------------------------------------------------- |
---|
233 | * QRsol solves the linear least squares problem |
---|
234 | * |
---|
235 | * min (b - H*x, b - H*x), x in R^n, |
---|
236 | * |
---|
237 | * where H is a Hessenberg matrix, and b is in R^(n+1). |
---|
238 | * It uses the QR factors of H computed by QRfact. |
---|
239 | * |
---|
240 | * n is the problem size; the matrix H is (n+1) by n. |
---|
241 | * |
---|
242 | * h is a matrix (computed by QRfact) containing the upper |
---|
243 | * triangular factor R of the original Hessenberg matrix H. |
---|
244 | * |
---|
245 | * q is an array of length 2*n (computed by QRfact) containing |
---|
246 | * the Givens rotations used to factor H. |
---|
247 | * |
---|
248 | * b is the (n+1)-vector appearing in the least squares problem |
---|
249 | * above. |
---|
250 | * |
---|
251 | * On return, b contains the solution x of the least squares |
---|
252 | * problem, if QRsol was successful. |
---|
253 | * |
---|
254 | * QRsol returns a 0 if successful. Otherwise, a zero was |
---|
255 | * encountered on the diagonal of the triangular factor R. |
---|
256 | * In this case, QRsol returns the equation number (numbered |
---|
257 | * from 1, not 0) of the zero entry. |
---|
258 | * ----------------------------------------------------------------- |
---|
259 | */ |
---|
260 | |
---|
261 | SUNDIALS_EXPORT int QRsol(int n, realtype **h, realtype *q, realtype *b); |
---|
262 | |
---|
263 | #ifdef __cplusplus |
---|
264 | } |
---|
265 | #endif |
---|
266 | |
---|
267 | #endif |
---|