1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 4272 $
4 * $Date: 2014-12-02 11:19:41 -0800 (Tue, 02 Dec 2014) $
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 KINBAND linear solver.
19 * -----------------------------------------------------------------
20 */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include <kinsol/kinsol_band.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 /* KINBAND linit, lsetup, lsolve, and lfree routines */
44 static int kinBandInit(KINMem kin_mem);
45 static int kinBandSetup(KINMem kin_mem);
46 static int kinBandsolve(KINMem kin_mem, N_Vector x, N_Vector b,
47 realtype *sJpnorm, realtype *sFdotJp);
48 static void kinBandFree(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 bjac (kindls_mem->d_bjac)
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 * KINBand
102 * -----------------------------------------------------------------
103 * This routine initializes the memory record and sets various function
104 * fields specific to the band linear solver module. KINBand first calls
105 * the existing lfree routine if this is not NULL. It then sets the
106 * cv_linit, cv_lsetup, cv_lsolve, and cv_lfree fields in (*cvode_mem)
107 * to be kinBandInit, kinBandSetup, kinBandsolve, and kinBandFree,
108 * respectively. It allocates memory for a structure of type
109 * KINDlsMemRec and sets the cv_lmem field in (*cvode_mem) to the
110 * address of this structure. It sets setupNonNull in (*cvode_mem) to be
111 * TRUE, b_mu to be mupper, b_ml to be mlower, and the bjac field to be
112 * kinDlsBandDQJac.
113 * Finally, it allocates memory for M, savedJ, and pivot.
114 *
115 * NOTE: The band linear solver assumes a serial implementation
116 * of the NVECTOR package. Therefore, KINBand will first
117 * test for compatible a compatible N_Vector internal
118 * representation by checking that the function
119 * N_VGetArrayPointer exists.
120 * -----------------------------------------------------------------
121 */
122
KINBand(void * kinmem,long int N,long int mupper,long int mlower)123 int KINBand(void *kinmem, long int N, long int mupper, long int mlower)
124 {
125 KINMem kin_mem;
126 KINDlsMem kindls_mem;
127
128 /* Return immediately if kinmem is NULL */
129 if (kinmem == NULL) {
130 KINProcessError(NULL, KINDLS_MEM_NULL, "KINBAND", "KINBand", MSGD_KINMEM_NULL);
131 return(KINDLS_MEM_NULL);
132 }
133 kin_mem = (KINMem) kinmem;
134
135 /* Test if the NVECTOR package is compatible with the BAND solver */
136 if (vec_tmpl->ops->nvgetarraypointer == NULL) {
137 KINProcessError(kin_mem, KINDLS_ILL_INPUT, "KINBAND", "KINBand", 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 = kinBandInit;
145 lsetup = kinBandSetup;
146 lsolve = kinBandsolve;
147 lfree = kinBandFree;
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, "KINBAND", "KINBand", MSGD_MEM_FAIL);
154 return(KINDLS_MEM_FAIL);
155 }
156
157 /* Set matrix type */
158 mtype = SUNDIALS_BAND;
159
160 /* Set default Jacobian routine and Jacobian data */
161 jacDQ = TRUE;
162 bjac = NULL;
163 J_data = NULL;
164 last_flag = KINDLS_SUCCESS;
165
166 setupNonNull = TRUE;
167
168 /* Load problem dimension */
169 n = N;
170
171 /* Load half-bandwiths in kindls_mem */
172 ml = mlower;
173 mu = mupper;
174
175 /* Test ml and mu for legality */
176 if ((ml < 0) || (mu < 0) || (ml >= N) || (mu >= N)) {
177 KINProcessError(kin_mem, KINDLS_MEM_FAIL, "KINBAND", "KINBand", MSGD_MEM_FAIL);
178 free(kindls_mem); kindls_mem = NULL;
179 return(KINDLS_ILL_INPUT);
180 }
181
182 /* Set extended upper half-bandwith for M (required for pivoting) */
183 smu = SUNMIN(N-1, mu + ml);
184
185 /* Allocate memory for J and pivot array */
186 J = NULL;
187 J = NewBandMat(N, mu, ml, smu);
188 if (J == NULL) {
189 KINProcessError(kin_mem, KINDLS_MEM_FAIL, "KINBAND", "KINBand", MSGD_MEM_FAIL);
190 free(kindls_mem); kindls_mem = NULL;
191 return(KINDLS_MEM_FAIL);
192 }
193
194 lpivots = NULL;
195 lpivots = NewLintArray(N);
196 if (lpivots == NULL) {
197 KINProcessError(kin_mem, KINDLS_MEM_FAIL, "KINBAND", "KINBand", MSGD_MEM_FAIL);
198 DestroyMat(J);
199 free(kindls_mem); kindls_mem = NULL;
200 return(KINDLS_MEM_FAIL);
201 }
202
203 /* This is a direct linear solver */
204 inexact_ls = FALSE;
205
206 /* Attach linear solver memory to integrator memory */
207 lmem = kindls_mem;
208
209 return(KINDLS_SUCCESS);
210 }
211
212 /*
213 * =================================================================
214 * PRIVATE FUNCTIONS
215 * =================================================================
216 */
217
218 /*
219 * -----------------------------------------------------------------
220 * kinBandInit
221 * -----------------------------------------------------------------
222 * This routine does remaining initializations specific to the band
223 * linear solver.
224 * -----------------------------------------------------------------
225 */
226
kinBandInit(KINMem kin_mem)227 static int kinBandInit(KINMem kin_mem)
228 {
229 KINDlsMem kindls_mem;
230
231 kindls_mem = (KINDlsMem) lmem;
232
233 nje = 0;
234 nfeDQ = 0;
235
236 if (jacDQ) {
237 bjac = kinDlsBandDQJac;
238 J_data = kin_mem;
239 } else {
240 J_data = kin_mem->kin_user_data;
241 }
242
243 /* Stop Picard solve if user fails to provide Jacobian */
244 if ( (strategy == KIN_PICARD) && jacDQ )
245 {
246 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "kinBandInit",
247 MSG_NOL_FAIL);
248 return(KIN_ILL_INPUT);
249 }
250
251 last_flag = KINDLS_SUCCESS;
252 return(0);
253 }
254
255 /*
256 * -----------------------------------------------------------------
257 * kinBandSetup
258 * -----------------------------------------------------------------
259 * This routine does the setup operations for the band linear solver.
260 * It makes a decision whether or not to call the Jacobian evaluation
261 * routine based on various state variables, and if not it uses the
262 * saved copy. In any case, it constructs the Newton matrix J,
263 * updates counters, and calls the band LU factorization routine.
264 * -----------------------------------------------------------------
265 */
266
kinBandSetup(KINMem kin_mem)267 static int kinBandSetup(KINMem kin_mem)
268 {
269 KINDlsMem kindls_mem;
270 int retval;
271 long int ier;
272
273 kindls_mem = (KINDlsMem) lmem;
274
275 nje++;
276 SetToZero(J);
277 retval = bjac(n, mu, ml, uu, fval, J, J_data, vtemp1, vtemp2);
278 if (retval != 0) {
279 last_flag = -1;
280 return(-1);
281 }
282
283 /* Do LU factorization of J */
284 ier = BandGBTRF(J, lpivots);
285
286 /* Return 0 if the LU was complete; otherwise return -1 */
287 last_flag = ier;
288 if (ier > 0) return(-1);
289
290 return(0);
291 }
292
293 /*
294 * -----------------------------------------------------------------
295 * kinBandsolve
296 * -----------------------------------------------------------------
297 * This routine handles the solve operation for the band linear solver
298 * by calling the band backsolve routine. The return value is 0.
299 * The argument *sJpnorm is ignored.
300 * -----------------------------------------------------------------
301 */
302
kinBandsolve(KINMem kin_mem,N_Vector x,N_Vector b,realtype * sJpnorm,realtype * sFdotJp)303 static int kinBandsolve(KINMem kin_mem, N_Vector x, N_Vector b,
304 realtype *sJpnorm, realtype *sFdotJp)
305 {
306 KINDlsMem kindls_mem;
307 realtype *xd;
308
309 kindls_mem = (KINDlsMem) lmem;
310
311 /* Copy the right-hand side into x */
312
313 N_VScale(ONE, b, x);
314
315 xd = N_VGetArrayPointer(x);
316
317 /* Back-solve and get solution in x */
318
319 BandGBTRS(J, lpivots, xd);
320
321 /* Compute the term sFdotJp for use in the linesearch routine.
322 This term is subsequently corrected if the step is reduced by
323 constraints or the linesearch.
324
325 sFdotJp is the dot product of the scaled f vector and the scaled
326 vector J*p, where the scaling uses fscale. */
327
328 N_VProd(b, fscale, b);
329 N_VProd(b, fscale, b);
330 *sFdotJp = N_VDotProd(fval, b);
331
332 last_flag = KINDLS_SUCCESS;
333 return(0);
334 }
335
336 /*
337 * -----------------------------------------------------------------
338 * kinBandFree
339 * -----------------------------------------------------------------
340 * This routine frees memory specific to the band linear solver.
341 * -----------------------------------------------------------------
342 */
343
kinBandFree(KINMem kin_mem)344 static void kinBandFree(KINMem kin_mem)
345 {
346 KINDlsMem kindls_mem;
347
348 kindls_mem = (KINDlsMem) lmem;
349
350 DestroyMat(J);
351 DestroyArray(lpivots);
352 free(kindls_mem); kindls_mem = NULL;
353 }
354
355