1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 1.11 $
4 * $Date: 2010/12/01 22:43:33 $
5 * -----------------------------------------------------------------
6 * Programmer(s): Radu Serban @ LLNL
7 * -----------------------------------------------------------------
8 * Copyright (c) 2002, The Regents of the University of California.
9 * Produced at the Lawrence Livermore National Laboratory.
10 * All rights reserved.
11 * For details, see the LICENSE file.
12 * -----------------------------------------------------------------
13 * This is the implementation file for the KINDENSE linear solver.
14 * -----------------------------------------------------------------
15 */
16
17 #include <stdio.h>
18 #include <stdlib.h>
19
20 #include <kinsol/kinsol_dense.h>
21 #include "kinsol_direct_impl.h"
22 #include "kinsol_impl.h"
23
24 #include <sundials/sundials_math.h>
25
26 /* Constants */
27
28 #define ZERO RCONST(0.0)
29 #define ONE RCONST(1.0)
30 #define TWO RCONST(2.0)
31
32 /*
33 * =================================================================
34 * PROTOTYPES FOR PRIVATE FUNCTIONS
35 * =================================================================
36 */
37
38 /* KINDENSE linit, lsetup, lsolve, and lfree routines */
39 static int kinDenseInit(KINMem kin_mem);
40 static int kinDenseSetup(KINMem kin_mem);
41 static int kinDenseSolve(KINMem kin_mem, N_Vector x, N_Vector b,
42 realtype *res_norm);
43 static void kinDenseFree(KINMem kin_mem);
44
45 /*
46 * =================================================================
47 * READIBILITY REPLACEMENTS
48 * =================================================================
49 */
50
51 #define lrw1 (kin_mem->kin_lrw1)
52 #define liw1 (kin_mem->kin_liw1)
53 #define func (kin_mem->kin_func)
54 #define printfl (kin_mem->kin_printfl)
55 #define linit (kin_mem->kin_linit)
56 #define lsetup (kin_mem->kin_lsetup)
57 #define lsolve (kin_mem->kin_lsolve)
58 #define lfree (kin_mem->kin_lfree)
59 #define lmem (kin_mem->kin_lmem)
60 #define inexact_ls (kin_mem->kin_inexact_ls)
61 #define uu (kin_mem->kin_uu)
62 #define fval (kin_mem->kin_fval)
63 #define uscale (kin_mem->kin_uscale)
64 #define fscale (kin_mem->kin_fscale)
65 #define sqrt_relfunc (kin_mem->kin_sqrt_relfunc)
66 #define sJpnorm (kin_mem->kin_sJpnorm)
67 #define sfdotJp (kin_mem->kin_sfdotJp)
68 #define errfp (kin_mem->kin_errfp)
69 #define infofp (kin_mem->kin_infofp)
70 #define setupNonNull (kin_mem->kin_setupNonNull)
71 #define vtemp1 (kin_mem->kin_vtemp1)
72 #define vec_tmpl (kin_mem->kin_vtemp1)
73 #define vtemp2 (kin_mem->kin_vtemp2)
74
75 #define mtype (kindls_mem->d_type)
76 #define n (kindls_mem->d_n)
77 #define ml (kindls_mem->d_ml)
78 #define mu (kindls_mem->d_mu)
79 #define smu (kindls_mem->d_smu)
80 #define jacDQ (kindls_mem->d_jacDQ)
81 #define djac (kindls_mem->d_djac)
82 #define J (kindls_mem->d_J)
83 #define lpivots (kindls_mem->d_lpivots)
84 #define nje (kindls_mem->d_nje)
85 #define nfeDQ (kindls_mem->d_nfeDQ)
86 #define J_data (kindls_mem->d_J_data)
87 #define last_flag (kindls_mem->d_last_flag)
88
89 /*
90 * =================================================================
91 * EXPORTED FUNCTIONS
92 * =================================================================
93 */
94
95 /*
96 * -----------------------------------------------------------------
97 * KINDense
98 * -----------------------------------------------------------------
99 * This routine initializes the memory record and sets various function
100 * fields specific to the dense linear solver module.
101 * KINDense sets the kin_linit, kin_lsetup, kin_lsolve, kin_lfree fields
102 * in *kinmem to be kinDenseInit, kinDenseSetup, kinDenseSolve, and
103 * kinDenseFree, respectively.
104 * It allocates memory for a structure of type KINDlsMemRec and sets
105 * the kin_lmem field in *kinmem to the address of this structure.
106 * It sets setupNonNull in *kinmem to TRUE, and the djac field to the
107 * default kinDlsDenseDQJac.
108 * Finally, it allocates memory for J and lpivots.
109 *
110 * NOTE: The dense linear solver assumes a serial implementation
111 * of the NVECTOR package. Therefore, KINDense will first
112 * test for compatible a compatible N_Vector internal
113 * representation by checking that N_VGetArrayPointer and
114 * N_VSetArrayPointer exist.
115 * -----------------------------------------------------------------
116 */
117
KINDense(void * kinmem,long int N)118 int KINDense(void *kinmem, long int N)
119 {
120 KINMem kin_mem;
121 KINDlsMem kindls_mem;
122
123 /* Return immediately if kinmem is NULL */
124 if (kinmem == NULL)
125 {
126 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDENSE", "KINDense", MSGD_KINMEM_NULL);
127 return(KINDLS_MEM_NULL);
128 }
129 kin_mem = (KINMem) kinmem;
130
131 /* Test if the NVECTOR package is compatible with the DENSE solver */
132 if (vec_tmpl->ops->nvgetarraypointer == NULL ||
133 vec_tmpl->ops->nvsetarraypointer == NULL)
134 {
135 KINProcessError(kin_mem, KINDLS_ILL_INPUT, "KINDENSE", "KINDense", MSGD_BAD_NVECTOR);
136 return(KINDLS_ILL_INPUT);
137 }
138
139 if (lfree != NULL)
140 {
141 lfree(kin_mem);
142 }
143
144 /* Set four main function fields in kin_mem */
145 linit = kinDenseInit;
146 lsetup = kinDenseSetup;
147 lsolve = kinDenseSolve;
148 lfree = kinDenseFree;
149
150 /* Get memory for KINDlsMemRec */
151 kindls_mem = NULL;
152 kindls_mem = (KINDlsMem) malloc(sizeof(struct KINDlsMemRec));
153 if (kindls_mem == NULL)
154 {
155 KINProcessError(kin_mem, KINDLS_MEM_FAIL, "KINDENSE", "KINDense", MSGD_MEM_FAIL);
156 return(KINDLS_MEM_FAIL);
157 }
158
159 /* Set matrix type */
160 mtype = SUNDIALS_DENSE;
161
162 /* Set default Jacobian routine and Jacobian data */
163 jacDQ = TRUE;
164 djac = NULL;
165 J_data = NULL;
166 last_flag = KINDLS_SUCCESS;
167
168 setupNonNull = TRUE;
169
170 /* Set problem dimension */
171 n = N;
172
173 /* Allocate memory for J and pivot array */
174
175 J = NULL;
176 J = NewDenseMat(N, N);
177 if (J == NULL)
178 {
179 KINProcessError(kin_mem, KINDLS_MEM_FAIL, "KINDENSE", "KINDense", MSGD_MEM_FAIL);
180 free(kindls_mem);
181 kindls_mem = NULL;
182 return(KINDLS_MEM_FAIL);
183 }
184
185 lpivots = NULL;
186 lpivots = NewLintArray(N);
187 if (lpivots == NULL)
188 {
189 KINProcessError(kin_mem, KINDLS_MEM_FAIL, "KINDENSE", "KINDense", MSGD_MEM_FAIL);
190 DestroyMat(J);
191 free(kindls_mem);
192 kindls_mem = NULL;
193 return(KINDLS_MEM_FAIL);
194 }
195
196 /* This is a direct linear solver */
197 inexact_ls = FALSE;
198
199 /* Attach linear solver memory to integrator memory */
200 lmem = kindls_mem;
201
202 return(KINDLS_SUCCESS);
203 }
204
205 /*
206 * =================================================================
207 * PRIVATE FUNCTIONS
208 * =================================================================
209 */
210
211 /*
212 * -----------------------------------------------------------------
213 * kinDenseInit
214 * -----------------------------------------------------------------
215 * This routine does remaining initializations specific to the dense
216 * linear solver.
217 * -----------------------------------------------------------------
218 */
219
kinDenseInit(KINMem kin_mem)220 static int kinDenseInit(KINMem kin_mem)
221 {
222 KINDlsMem kindls_mem;
223
224 kindls_mem = (KINDlsMem) lmem;
225
226 nje = 0;
227 nfeDQ = 0;
228
229 if (jacDQ)
230 {
231 djac = kinDlsDenseDQJac;
232 J_data = kin_mem;
233 }
234 else
235 {
236 J_data = kin_mem->kin_user_data;
237 }
238
239 last_flag = KINDLS_SUCCESS;
240 return(0);
241 }
242
243 /*
244 * -----------------------------------------------------------------
245 * kinDenseSetup
246 * -----------------------------------------------------------------
247 * This routine does the setup operations for the dense linear solver.
248 * It calls the dense LU factorization routine.
249 * -----------------------------------------------------------------
250 */
251
kinDenseSetup(KINMem kin_mem)252 static int kinDenseSetup(KINMem kin_mem)
253 {
254 KINDlsMem kindls_mem;
255 int retval;
256 long int ier;
257
258 kindls_mem = (KINDlsMem) lmem;
259
260 nje++;
261 SetToZero(J);
262 retval = djac(n, uu, fval, J, J_data, vtemp1, vtemp2);
263 if (retval != 0)
264 {
265 last_flag = -1;
266 return(-1);
267 }
268
269 /* Do LU factorization of J */
270 ier = DenseGETRF(J, lpivots);
271
272 /* Return 0 if the LU was complete; otherwise return -1 */
273 last_flag = ier;
274 if (ier > 0)
275 {
276 return(-1);
277 }
278
279 return(0);
280 }
281
282 /*
283 * -----------------------------------------------------------------
284 * kinDenseSolve
285 * -----------------------------------------------------------------
286 * This routine handles the solve operation for the dense linear solver
287 * by calling the dense backsolve routine. The returned value is 0.
288 * -----------------------------------------------------------------
289 */
290
kinDenseSolve(KINMem kin_mem,N_Vector x,N_Vector b,realtype * res_norm)291 static int kinDenseSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm)
292 {
293 KINDlsMem kindls_mem;
294 realtype *xd;
295
296 kindls_mem = (KINDlsMem) lmem;
297
298 /* Copy the right-hand side into x */
299
300 N_VScale(ONE, b, x);
301
302 xd = N_VGetArrayPointer(x);
303
304 /* Back-solve and get solution in x */
305
306 DenseGETRS(J, lpivots, xd);
307
308 /* Compute the terms Jpnorm and sfdotJp for use in the global strategy
309 routines and in KINForcingTerm. Both of these terms are subsequently
310 corrected if the step is reduced by constraints or the line search.
311
312 sJpnorm is the norm of the scaled product (scaled by fscale) of
313 the current Jacobian matrix J and the step vector p.
314
315 sfdotJp is the dot product of the scaled f vector and the scaled
316 vector J*p, where the scaling uses fscale. */
317
318 sJpnorm = N_VWL2Norm(b, fscale);
319 N_VProd(b, fscale, b);
320 N_VProd(b, fscale, b);
321 sfdotJp = N_VDotProd(fval, b);
322
323 last_flag = KINDLS_SUCCESS;
324
325 return(0);
326 }
327
328 /*
329 * -----------------------------------------------------------------
330 * kinDenseFree
331 * -----------------------------------------------------------------
332 * This routine frees memory specific to the dense linear solver.
333 * -----------------------------------------------------------------
334 */
335
kinDenseFree(KINMem kin_mem)336 static void kinDenseFree(KINMem kin_mem)
337 {
338 KINDlsMem kindls_mem;
339
340 kindls_mem = (KINDlsMem) lmem;
341
342 DestroyMat(J);
343 DestroyArray(lpivots);
344 free(kindls_mem);
345 kindls_mem = NULL;
346 }
347
348