1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 1.5 $
4 * $Date: 2010/12/01 22:21:04 $
5 * -----------------------------------------------------------------
6 * Programmer: Radu Serban @ LLNL
7 * -----------------------------------------------------------------
8 * Copyright (c) 2006, 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 CVDLS linear solvers
14 * -----------------------------------------------------------------
15 */
16
17 /*
18 * =================================================================
19 * IMPORTED HEADER FILES
20 * =================================================================
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "cvode_impl.h"
27 #include "cvode_direct_impl.h"
28 #include <sundials/sundials_math.h>
29
30 /*
31 * =================================================================
32 * FUNCTION SPECIFIC CONSTANTS
33 * =================================================================
34 */
35
36 /* Constant for DQ Jacobian approximation */
37 #define MIN_INC_MULT RCONST(1000.0)
38
39 #define ZERO RCONST(0.0)
40 #define ONE RCONST(1.0)
41 #define TWO RCONST(2.0)
42
43 /*
44 * =================================================================
45 * READIBILITY REPLACEMENTS
46 * =================================================================
47 */
48
49 #define f (cv_mem->cv_f)
50 #define user_data (cv_mem->cv_user_data)
51 #define uround (cv_mem->cv_uround)
52 #define nst (cv_mem->cv_nst)
53 #define tn (cv_mem->cv_tn)
54 #define h (cv_mem->cv_h)
55 #define gamma (cv_mem->cv_gamma)
56 #define gammap (cv_mem->cv_gammap)
57 #define gamrat (cv_mem->cv_gamrat)
58 #define ewt (cv_mem->cv_ewt)
59
60 #define lmem (cv_mem->cv_lmem)
61
62 #define mtype (cvdls_mem->d_type)
63 #define n (cvdls_mem->d_n)
64 #define ml (cvdls_mem->d_ml)
65 #define mu (cvdls_mem->d_mu)
66 #define smu (cvdls_mem->d_smu)
67 #define jacDQ (cvdls_mem->d_jacDQ)
68 #define djac (cvdls_mem->d_djac)
69 #define bjac (cvdls_mem->d_bjac)
70 #define M (cvdls_mem->d_M)
71 #define nje (cvdls_mem->d_nje)
72 #define nfeDQ (cvdls_mem->d_nfeDQ)
73 #define last_flag (cvdls_mem->d_last_flag)
74
75 /*
76 * =================================================================
77 * EXPORTED FUNCTIONS
78 * =================================================================
79 */
80
81 /*
82 * CVDlsSetDenseJacFn specifies the dense Jacobian function.
83 */
CVDlsSetDenseJacFn(void * cvode_mem,CVDlsDenseJacFn jac)84 int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn jac)
85 {
86 CVodeMem cv_mem;
87 CVDlsMem cvdls_mem;
88
89 /* Return immediately if cvode_mem is NULL */
90 if (cvode_mem == NULL)
91 {
92 CVProcessError(NULL, CVDLS_MEM_NULL, "CVDLS", "CVDlsSetDenseJacFn", MSGD_CVMEM_NULL);
93 return(CVDLS_MEM_NULL);
94 }
95 cv_mem = (CVodeMem) cvode_mem;
96
97 if (lmem == NULL)
98 {
99 CVProcessError(cv_mem, CVDLS_LMEM_NULL, "CVDLS", "CVDlsSetDenseJacFn", MSGD_LMEM_NULL);
100 return(CVDLS_LMEM_NULL);
101 }
102 cvdls_mem = (CVDlsMem) lmem;
103
104 if (jac != NULL)
105 {
106 jacDQ = FALSE;
107 djac = jac;
108 }
109 else
110 {
111 jacDQ = TRUE;
112 }
113
114 return(CVDLS_SUCCESS);
115 }
116
117 /*
118 * CVDlsSetBandJacFn specifies the band Jacobian function.
119 */
CVDlsSetBandJacFn(void * cvode_mem,CVDlsBandJacFn jac)120 int CVDlsSetBandJacFn(void *cvode_mem, CVDlsBandJacFn jac)
121 {
122 CVodeMem cv_mem;
123 CVDlsMem cvdls_mem;
124
125 /* Return immediately if cvode_mem is NULL */
126 if (cvode_mem == NULL)
127 {
128 CVProcessError(NULL, CVDLS_MEM_NULL, "CVDLS", "CVDlsSetBandJacFn", MSGD_CVMEM_NULL);
129 return(CVDLS_MEM_NULL);
130 }
131 cv_mem = (CVodeMem) cvode_mem;
132
133 if (lmem == NULL)
134 {
135 CVProcessError(cv_mem, CVDLS_LMEM_NULL, "CVDLS", "CVDlsSetBandJacFn", MSGD_LMEM_NULL);
136 return(CVDLS_LMEM_NULL);
137 }
138 cvdls_mem = (CVDlsMem) lmem;
139
140 if (jac != NULL)
141 {
142 jacDQ = FALSE;
143 bjac = jac;
144 }
145 else
146 {
147 jacDQ = TRUE;
148 }
149
150 return(CVDLS_SUCCESS);
151 }
152
153 /*
154 * CVDlsGetWorkSpace returns the length of workspace allocated for the
155 * CVDLS linear solver.
156 */
CVDlsGetWorkSpace(void * cvode_mem,long int * lenrwLS,long int * leniwLS)157 int CVDlsGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS)
158 {
159 CVodeMem cv_mem;
160 CVDlsMem cvdls_mem;
161
162 /* Return immediately if cvode_mem is NULL */
163 if (cvode_mem == NULL)
164 {
165 CVProcessError(NULL, CVDLS_MEM_NULL, "CVDLS", "CVDlsGetWorkSpace", MSGD_CVMEM_NULL);
166 return(CVDLS_MEM_NULL);
167 }
168 cv_mem = (CVodeMem) cvode_mem;
169
170 if (lmem == NULL)
171 {
172 CVProcessError(cv_mem, CVDLS_LMEM_NULL, "CVDLS", "CVDlsGetWorkSpace", MSGD_LMEM_NULL);
173 return(CVDLS_LMEM_NULL);
174 }
175 cvdls_mem = (CVDlsMem) lmem;
176
177 if (mtype == SUNDIALS_DENSE)
178 {
179 *lenrwLS = 2 * n * n;
180 *leniwLS = n;
181 }
182 else if (mtype == SUNDIALS_BAND)
183 {
184 *lenrwLS = n * (smu + mu + 2 * ml + 2);
185 *leniwLS = n;
186 }
187
188 return(CVDLS_SUCCESS);
189 }
190
191 /*
192 * CVDlsGetNumJacEvals returns the number of Jacobian evaluations.
193 */
CVDlsGetNumJacEvals(void * cvode_mem,long int * njevals)194 int CVDlsGetNumJacEvals(void *cvode_mem, long int *njevals)
195 {
196 CVodeMem cv_mem;
197 CVDlsMem cvdls_mem;
198
199 /* Return immediately if cvode_mem is NULL */
200 if (cvode_mem == NULL)
201 {
202 CVProcessError(NULL, CVDLS_MEM_NULL, "CVDLS", "CVDlsGetNumJacEvals", MSGD_CVMEM_NULL);
203 return(CVDLS_MEM_NULL);
204 }
205 cv_mem = (CVodeMem) cvode_mem;
206
207 if (lmem == NULL)
208 {
209 CVProcessError(cv_mem, CVDLS_LMEM_NULL, "CVDLS", "CVDlsGetNumJacEvals", MSGD_LMEM_NULL);
210 return(CVDLS_LMEM_NULL);
211 }
212 cvdls_mem = (CVDlsMem) lmem;
213
214 *njevals = nje;
215
216 return(CVDLS_SUCCESS);
217 }
218
219 /*
220 * CVDlsGetNumRhsEvals returns the number of calls to the ODE function
221 * needed for the DQ Jacobian approximation.
222 */
CVDlsGetNumRhsEvals(void * cvode_mem,long int * nfevalsLS)223 int CVDlsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS)
224 {
225 CVodeMem cv_mem;
226 CVDlsMem cvdls_mem;
227
228 /* Return immediately if cvode_mem is NULL */
229 if (cvode_mem == NULL)
230 {
231 CVProcessError(NULL, CVDLS_MEM_NULL, "CVDLS", "CVDlsGetNumRhsEvals", MSGD_CVMEM_NULL);
232 return(CVDLS_MEM_NULL);
233 }
234 cv_mem = (CVodeMem) cvode_mem;
235
236 if (lmem == NULL)
237 {
238 CVProcessError(cv_mem, CVDLS_LMEM_NULL, "CVDLS", "CVDlsGetNumRhsEvals", MSGD_LMEM_NULL);
239 return(CVDLS_LMEM_NULL);
240 }
241 cvdls_mem = (CVDlsMem) lmem;
242
243 *nfevalsLS = nfeDQ;
244
245 return(CVDLS_SUCCESS);
246 }
247
248 /*
249 * CVDlsGetReturnFlagName returns the name associated with a CVDLS
250 * return value.
251 */
CVDlsGetReturnFlagName(long int flag)252 char *CVDlsGetReturnFlagName(long int flag)
253 {
254 char *name;
255
256 name = (char *)malloc(30 * sizeof(char));
257
258 switch (flag)
259 {
260 case CVDLS_SUCCESS:
261 sprintf(name, "CVDLS_SUCCESS");
262 break;
263 case CVDLS_MEM_NULL:
264 sprintf(name, "CVDLS_MEM_NULL");
265 break;
266 case CVDLS_LMEM_NULL:
267 sprintf(name, "CVDLS_LMEM_NULL");
268 break;
269 case CVDLS_ILL_INPUT:
270 sprintf(name, "CVDLS_ILL_INPUT");
271 break;
272 case CVDLS_MEM_FAIL:
273 sprintf(name, "CVDLS_MEM_FAIL");
274 break;
275 case CVDLS_JACFUNC_UNRECVR:
276 sprintf(name, "CVDLS_JACFUNC_UNRECVR");
277 break;
278 case CVDLS_JACFUNC_RECVR:
279 sprintf(name, "CVDLS_JACFUNC_RECVR");
280 break;
281 default:
282 sprintf(name, "NONE");
283 }
284
285 return(name);
286 }
287
288 /*
289 * CVDlsGetLastFlag returns the last flag set in a CVDLS function.
290 */
CVDlsGetLastFlag(void * cvode_mem,long int * flag)291 int CVDlsGetLastFlag(void *cvode_mem, long int *flag)
292 {
293 CVodeMem cv_mem;
294 CVDlsMem cvdls_mem;
295
296 /* Return immediately if cvode_mem is NULL */
297 if (cvode_mem == NULL)
298 {
299 CVProcessError(NULL, CVDLS_MEM_NULL, "CVDLS", "CVDlsGetLastFlag", MSGD_CVMEM_NULL);
300 return(CVDLS_MEM_NULL);
301 }
302 cv_mem = (CVodeMem) cvode_mem;
303
304 if (lmem == NULL)
305 {
306 CVProcessError(cv_mem, CVDLS_LMEM_NULL, "CVDLS", "CVDlsGetLastFlag", MSGD_LMEM_NULL);
307 return(CVDLS_LMEM_NULL);
308 }
309 cvdls_mem = (CVDlsMem) lmem;
310
311 *flag = last_flag;
312
313 return(CVDLS_SUCCESS);
314 }
315
316 /*
317 * =================================================================
318 * DQ JACOBIAN APPROXIMATIONS
319 * =================================================================
320 */
321
322 /*
323 * -----------------------------------------------------------------
324 * cvDlsDenseDQJac
325 * -----------------------------------------------------------------
326 * This routine generates a dense difference quotient approximation to
327 * the Jacobian of f(t,y). It assumes that a dense matrix of type
328 * DlsMat is stored column-wise, and that elements within each column
329 * are contiguous. The address of the jth column of J is obtained via
330 * the macro DENSE_COL and this pointer is associated with an N_Vector
331 * using the N_VGetArrayPointer/N_VSetArrayPointer functions.
332 * Finally, the actual computation of the jth column of the Jacobian is
333 * done with a call to N_VLinearSum.
334 * -----------------------------------------------------------------
335 */
336
cvDlsDenseDQJac(long int N,realtype t,N_Vector y,N_Vector fy,DlsMat Jac,void * data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)337 int cvDlsDenseDQJac(long int N, realtype t,
338 N_Vector y, N_Vector fy,
339 DlsMat Jac, void *data,
340 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
341 {
342 realtype fnorm, minInc, inc, inc_inv, yjsaved, srur;
343 realtype *tmp2_data, *y_data, *ewt_data;
344 N_Vector ftemp, jthCol;
345 long int j;
346 int retval = 0;
347
348 CVodeMem cv_mem;
349 CVDlsMem cvdls_mem;
350
351 /* data points to cvode_mem */
352 cv_mem = (CVodeMem) data;
353 cvdls_mem = (CVDlsMem) lmem;
354
355 /* Save pointer to the array in tmp2 */
356 tmp2_data = N_VGetArrayPointer(tmp2);
357
358 /* Rename work vectors for readibility */
359 ftemp = tmp1;
360 jthCol = tmp2;
361
362 /* Obtain pointers to the data for ewt, y */
363 ewt_data = N_VGetArrayPointer(ewt);
364 y_data = N_VGetArrayPointer(y);
365
366 /* Set minimum increment based on uround and norm of f */
367 srur = RSqrt(uround);
368 fnorm = N_VWrmsNorm(fy, ewt);
369 minInc = (fnorm != ZERO) ?
370 (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
371
372 for (j = 0; j < N; j++)
373 {
374
375 /* Generate the jth col of J(tn,y) */
376
377 N_VSetArrayPointer(DENSE_COL(Jac, j), jthCol);
378
379 yjsaved = y_data[j];
380 inc = MAX(srur * ABS(yjsaved), minInc / ewt_data[j]);
381 y_data[j] += inc;
382
383 retval = f(t, y, ftemp, user_data);
384 nfeDQ++;
385 if (retval != 0)
386 {
387 break;
388 }
389
390 y_data[j] = yjsaved;
391
392 inc_inv = ONE / inc;
393 N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);
394
395 DENSE_COL(Jac, j) = N_VGetArrayPointer(jthCol);
396 }
397
398 /* Restore original array pointer in tmp2 */
399 N_VSetArrayPointer(tmp2_data, tmp2);
400
401 return(retval);
402 }
403
404 /*
405 * -----------------------------------------------------------------
406 * cvDlsBandDQJac
407 * -----------------------------------------------------------------
408 * This routine generates a banded difference quotient approximation to
409 * the Jacobian of f(t,y). It assumes that a band matrix of type
410 * DlsMat is stored column-wise, and that elements within each column
411 * are contiguous. This makes it possible to get the address of a column
412 * of J via the macro BAND_COL and to write a simple for loop to set
413 * each of the elements of a column in succession.
414 * -----------------------------------------------------------------
415 */
416
cvDlsBandDQJac(long int N,long int mupper,long int mlower,realtype t,N_Vector y,N_Vector fy,DlsMat Jac,void * data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)417 int cvDlsBandDQJac(long int N, long int mupper, long int mlower,
418 realtype t, N_Vector y, N_Vector fy,
419 DlsMat Jac, void *data,
420 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
421 {
422 N_Vector ftemp, ytemp;
423 realtype fnorm, minInc, inc, inc_inv, srur;
424 realtype *col_j, *ewt_data, *fy_data, *ftemp_data, *y_data, *ytemp_data;
425 long int group, i, j, width, ngroups, i1, i2;
426 int retval = 0;
427
428 CVodeMem cv_mem;
429 CVDlsMem cvdls_mem;
430
431 /* data points to cvode_mem */
432 cv_mem = (CVodeMem) data;
433 cvdls_mem = (CVDlsMem) lmem;
434
435 /* Rename work vectors for use as temporary values of y and f */
436 ftemp = tmp1;
437 ytemp = tmp2;
438
439 /* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp */
440 ewt_data = N_VGetArrayPointer(ewt);
441 fy_data = N_VGetArrayPointer(fy);
442 ftemp_data = N_VGetArrayPointer(ftemp);
443 y_data = N_VGetArrayPointer(y);
444 ytemp_data = N_VGetArrayPointer(ytemp);
445
446 /* Load ytemp with y = predicted y vector */
447 N_VScale(ONE, y, ytemp);
448
449 /* Set minimum increment based on uround and norm of f */
450 srur = RSqrt(uround);
451 fnorm = N_VWrmsNorm(fy, ewt);
452 minInc = (fnorm != ZERO) ?
453 (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
454
455 /* Set bandwidth and number of column groups for band differencing */
456 width = mlower + mupper + 1;
457 ngroups = MIN(width, N);
458
459 /* Loop over column groups. */
460 for (group = 1; group <= ngroups; group++)
461 {
462
463 /* Increment all y_j in group */
464 for (j = group - 1; j < N; j += width)
465 {
466 inc = MAX(srur * ABS(y_data[j]), minInc / ewt_data[j]);
467 ytemp_data[j] += inc;
468 }
469
470 /* Evaluate f with incremented y */
471
472 retval = f(tn, ytemp, ftemp, user_data);
473 nfeDQ++;
474 if (retval != 0)
475 {
476 break;
477 }
478
479 /* Restore ytemp, then form and load difference quotients */
480 for (j = group - 1; j < N; j += width)
481 {
482 ytemp_data[j] = y_data[j];
483 col_j = BAND_COL(Jac, j);
484 inc = MAX(srur * ABS(y_data[j]), minInc / ewt_data[j]);
485 inc_inv = ONE / inc;
486 i1 = MAX(0, j - mupper);
487 i2 = MIN(j + mlower, N - 1);
488 for (i = i1; i <= i2; i++)
489 {
490 BAND_COL_ELEM(col_j, i, j) = inc_inv * (ftemp_data[i] - fy_data[i]);
491 }
492 }
493 }
494
495 return(retval);
496 }
497
498