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