1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 1.1 $
4 * $Date: 2006/11/22 00:12:48 $
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 CPDIRECT 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 "cpodes_direct_impl.h"
27 #include "cpodes_private.h"
28
29 #include <sundials/sundials_math.h>
30
31 /*
32 * =================================================================
33 * FUNCTION SPECIFIC CONSTANTS
34 * =================================================================
35 */
36
37 /* Constant for DQ Jacobian approximation */
38 #define MIN_INC_MULT RCONST(1000.0)
39
40 #define ZERO RCONST(0.0)
41 #define ONE RCONST(1.0)
42 #define TWO RCONST(2.0)
43
44 /*
45 * =================================================================
46 * READIBILITY REPLACEMENTS
47 * =================================================================
48 */
49
50 #define ode_type (cp_mem->cp_ode_type)
51 #define fe (cp_mem->cp_fe)
52 #define fi (cp_mem->cp_fi)
53 #define f_data (cp_mem->cp_f_data)
54 #define uround (cp_mem->cp_uround)
55 #define nst (cp_mem->cp_nst)
56 #define tn (cp_mem->cp_tn)
57 #define h (cp_mem->cp_h)
58 #define gamma (cp_mem->cp_gamma)
59 #define gammap (cp_mem->cp_gammap)
60 #define gamrat (cp_mem->cp_gamrat)
61 #define ewt (cp_mem->cp_ewt)
62 #define lmem (cp_mem->cp_lmem)
63
64 #define cfun (cp_mem->cp_cfun)
65 #define c_data (cp_mem->cp_c_data)
66 #define lmemP (cp_mem->cp_lmemP)
67
68 #define mtype (cpdls_mem->d_type)
69 #define n (cpdls_mem->d_n)
70 #define ml (cpdls_mem->d_ml)
71 #define mu (cpdls_mem->d_mu)
72 #define smu (cpdls_mem->d_smu)
73 #define djacE (cpdls_mem->d_djacE)
74 #define djacI (cpdls_mem->d_djacI)
75 #define bjacE (cpdls_mem->d_bjacE)
76 #define bjacI (cpdls_mem->d_bjacI)
77 #define M (cpdls_mem->d_M)
78 #define savedJ (cpdls_mem->d_savedJ)
79 #define pivots (cpdls_mem->d_pivots)
80 #define nstlj (cpdls_mem->d_nstlj)
81 #define nje (cpdls_mem->d_nje)
82 #define nfeDQ (cpdls_mem->d_nfeDQ)
83 #define J_data (cpdls_mem->d_J_data)
84 #define last_flag (cpdls_mem->d_last_flag)
85
86 #define djacP (cpdlsP_mem->d_jacP)
87 #define JP_data (cpdlsP_mem->d_JP_data)
88 #define njeP (cpdlsP_mem->d_njeP)
89 #define nceDQ (cpdlsP_mem->d_nceDQ)
90
91 /*
92 * =================================================================
93 * EXPORTED FUNCTIONS FOR IMPLICIT INTEGRATION
94 * =================================================================
95 */
96
97 /*
98 * CPDlsSetJacFn specifies the (dense or band) Jacobian function.
99 */
CPDlsSetJacFn(void * cpode_mem,void * jac,void * jac_data)100 int CPDlsSetJacFn(void *cpode_mem, void *jac, void *jac_data)
101 {
102 CPodeMem cp_mem;
103 CPDlsMem cpdls_mem;
104
105 /* Return immediately if cpode_mem is NULL */
106 if (cpode_mem == NULL) {
107 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsSetJacFn", MSGD_CPMEM_NULL);
108 return(CPDIRECT_MEM_NULL);
109 }
110 cp_mem = (CPodeMem) cpode_mem;
111
112 if (lmem == NULL) {
113 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsSetJacFn", MSGD_LMEM_NULL);
114 return(CPDIRECT_LMEM_NULL);
115 }
116 cpdls_mem = (CPDlsMem) lmem;
117
118 switch (ode_type) {
119
120 case CP_EXPL:
121 if (mtype == SUNDIALS_DENSE)
122 djacE = (CPDlsDenseJacExplFn) jac;
123 else if (mtype == SUNDIALS_BAND)
124 bjacE = (CPDlsBandJacExplFn) jac;
125 break;
126
127 case CP_IMPL:
128 if (mtype == SUNDIALS_DENSE)
129 djacI = (CPDlsDenseJacImplFn) jac;
130 else if (mtype == SUNDIALS_BAND)
131 bjacI = (CPDlsBandJacImplFn) jac;
132 break;
133
134 }
135
136 J_data = jac_data;
137
138 return(CPDIRECT_SUCCESS);
139 }
140
141 /*
142 * CPDlsGetWorkSpace returns the length of workspace allocated for the
143 * CPDIRECT linear solver.
144 */
CPDlsGetWorkSpace(void * cpode_mem,long int * lenrwLS,long int * leniwLS)145 int CPDlsGetWorkSpace(void *cpode_mem, long int *lenrwLS, long int *leniwLS)
146 {
147 CPodeMem cp_mem;
148 CPDlsMem cpdls_mem;
149
150 /* Return immediately if cpode_mem is NULL */
151 if (cpode_mem == NULL) {
152 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsGetWorkSpace", MSGD_CPMEM_NULL);
153 return(CPDIRECT_MEM_NULL);
154 }
155 cp_mem = (CPodeMem) cpode_mem;
156
157 if (lmem == NULL) {
158 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsGetWorkSpace", MSGD_LMEM_NULL);
159 return(CPDIRECT_LMEM_NULL);
160 }
161 cpdls_mem = (CPDlsMem) lmem;
162
163 switch (ode_type) {
164
165 case CP_EXPL:
166
167 if (mtype == SUNDIALS_DENSE) {
168 *lenrwLS = 2*n*n;
169 *leniwLS = n;
170 } else if (mtype == SUNDIALS_BAND) {
171 *lenrwLS = n*(smu + mu + 2*ml + 2);
172 *leniwLS = n;
173 }
174
175 break;
176
177 case CP_IMPL:
178
179 if (mtype == SUNDIALS_DENSE) {
180 *lenrwLS = n*n;
181 *leniwLS = n;
182 } else if (mtype == SUNDIALS_BAND) {
183 *lenrwLS = n*(smu + ml + 2);
184 *leniwLS = n;
185 }
186
187 break;
188
189 }
190
191 return(CPDIRECT_SUCCESS);
192 }
193
194 /*
195 * CPDlsGetNumJacEvals returns the number of Jacobian evaluations.
196 */
CPDlsGetNumJacEvals(void * cpode_mem,long int * njevals)197 int CPDlsGetNumJacEvals(void *cpode_mem, long int *njevals)
198 {
199 CPodeMem cp_mem;
200 CPDlsMem cpdls_mem;
201
202 /* Return immediately if cpode_mem is NULL */
203 if (cpode_mem == NULL) {
204 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsGetNumJacEvals", MSGD_CPMEM_NULL);
205 return(CPDIRECT_MEM_NULL);
206 }
207 cp_mem = (CPodeMem) cpode_mem;
208
209 if (lmem == NULL) {
210 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsGetNumJacEvals", MSGD_LMEM_NULL);
211 return(CPDIRECT_LMEM_NULL);
212 }
213 cpdls_mem = (CPDlsMem) lmem;
214
215 *njevals = nje;
216
217 return(CPDIRECT_SUCCESS);
218 }
219
220 /*
221 * CPDlsGetNumFctEvals returns the number of calls to the ODE function
222 * needed for the DQ Jacobian approximation.
223 */
CPDlsGetNumFctEvals(void * cpode_mem,long int * nfevalsLS)224 int CPDlsGetNumFctEvals(void *cpode_mem, long int *nfevalsLS)
225 {
226 CPodeMem cp_mem;
227 CPDlsMem cpdls_mem;
228
229 /* Return immediately if cpode_mem is NULL */
230 if (cpode_mem == NULL) {
231 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsGetNumRhsEvals", MSGD_CPMEM_NULL);
232 return(CPDIRECT_MEM_NULL);
233 }
234 cp_mem = (CPodeMem) cpode_mem;
235
236 if (lmem == NULL) {
237 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsGetNumRhsEvals", MSGD_LMEM_NULL);
238 return(CPDIRECT_LMEM_NULL);
239 }
240 cpdls_mem = (CPDlsMem) lmem;
241
242 *nfevalsLS = nfeDQ;
243
244 return(CPDIRECT_SUCCESS);
245 }
246
247 /*
248 * CPDlsGetReturnFlagName returns the name associated with a CPDIRECT
249 * return value.
250 */
CPDlsGetReturnFlagName(int flag)251 char *CPDlsGetReturnFlagName(int flag)
252 {
253 char *name;
254
255 name = (char *)malloc(30*sizeof(char));
256
257 switch(flag) {
258 case CPDIRECT_SUCCESS:
259 sprintf(name,"CPDIRECT_SUCCESS");
260 break;
261 case CPDIRECT_MEM_NULL:
262 sprintf(name,"CPDIRECT_MEM_NULL");
263 break;
264 case CPDIRECT_LMEM_NULL:
265 sprintf(name,"CPDIRECT_LMEM_NULL");
266 break;
267 case CPDIRECT_ILL_INPUT:
268 sprintf(name,"CPDIRECT_ILL_INPUT");
269 break;
270 case CPDIRECT_MEM_FAIL:
271 sprintf(name,"CPDIRECT_MEM_FAIL");
272 break;
273 case CPDIRECT_JACFUNC_UNRECVR:
274 sprintf(name,"CPDIRECT_JACFUNC_UNRECVR");
275 break;
276 case CPDIRECT_JACFUNC_RECVR:
277 sprintf(name,"CPDIRECT_JACFUNC_RECVR");
278 break;
279 default:
280 sprintf(name,"NONE");
281 }
282
283 return(name);
284 }
285
286 /*
287 * CPDlsGetLastFlag returns the last flag set in a CPDIRECT function.
288 */
CPDlsGetLastFlag(void * cpode_mem,int * flag)289 int CPDlsGetLastFlag(void *cpode_mem, int *flag)
290 {
291 CPodeMem cp_mem;
292 CPDlsMem cpdls_mem;
293
294 /* Return immediately if cpode_mem is NULL */
295 if (cpode_mem == NULL) {
296 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsGetLastFlag", MSGD_CPMEM_NULL);
297 return(CPDIRECT_MEM_NULL);
298 }
299 cp_mem = (CPodeMem) cpode_mem;
300
301 if (lmem == NULL) {
302 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsGetLastFlag", MSGD_LMEM_NULL);
303 return(CPDIRECT_LMEM_NULL);
304 }
305 cpdls_mem = (CPDlsMem) lmem;
306
307 *flag = last_flag;
308
309 return(CPDIRECT_SUCCESS);
310 }
311
312 /*
313 * =================================================================
314 * EXPORTED FUNCTIONS FOR PROJECTION
315 * =================================================================
316 */
317
318 /*
319 * CPDlsProjSetJacFn specifies the constraint Jacobian function.
320 */
CPDlsProjSetJacFn(void * cpode_mem,void * jacP,void * jacP_data)321 int CPDlsProjSetJacFn(void *cpode_mem, void *jacP, void *jacP_data)
322 {
323 CPodeMem cp_mem;
324 CPDlsProjMem cpdlsP_mem;
325
326 /* Return immediately if cpode_mem is NULL */
327 if (cpode_mem == NULL) {
328 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsProjSetJacFn", MSGD_CPMEM_NULL);
329 return(CPDIRECT_MEM_NULL);
330 }
331 cp_mem = (CPodeMem) cpode_mem;
332
333 if (lmemP == NULL) {
334 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsProjSetJacFn", MSGD_LMEM_NULL);
335 return(CPDIRECT_LMEM_NULL);
336 }
337 cpdlsP_mem = (CPDlsProjMem) lmemP;
338
339 djacP = (CPDlsDenseProjJacFn) jacP;
340 JP_data = jacP_data;
341
342 return(CPDIRECT_SUCCESS);
343 }
344
345 /*
346 * CPDlsProjGetNumJacEvals returns the number of constraint Jacobian evaluations
347 */
CPDlsProjGetNumJacEvals(void * cpode_mem,long int * njPevals)348 int CPDlsProjGetNumJacEvals(void *cpode_mem, long int *njPevals)
349 {
350 CPodeMem cp_mem;
351 CPDlsProjMem cpdlsP_mem;
352
353 /* Return immediately if cpode_mem is NULL */
354 if (cpode_mem == NULL) {
355 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsProjGetNumJacEvals", MSGD_CPMEM_NULL);
356 return(CPDIRECT_MEM_NULL);
357 }
358 cp_mem = (CPodeMem) cpode_mem;
359
360 if (lmemP == NULL) {
361 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsProjGetNumJacEvals", MSGD_LMEM_NULL);
362 return(CPDIRECT_LMEM_NULL);
363 }
364 cpdlsP_mem = (CPDlsProjMem) lmemP;
365
366 *njPevals = njeP;
367
368 return(CPDIRECT_SUCCESS);
369 }
370
371 /*
372 * CPDlsProjGetNumFctEvals returns the number of constraint function
373 * evaluations for computing the DQ constraint Jacobian.
374 */
CPDlsProjGetNumFctEvals(void * cpode_mem,long int * ncevalsLS)375 int CPDlsProjGetNumFctEvals(void *cpode_mem, long int *ncevalsLS)
376 {
377 CPodeMem cp_mem;
378 CPDlsProjMem cpdlsP_mem;
379
380 /* Return immediately if cpode_mem is NULL */
381 if (cpode_mem == NULL) {
382 cpProcessError(NULL, CPDIRECT_MEM_NULL, "CPDIRECT", "CPDlsProjGetNumFctEvals", MSGD_CPMEM_NULL);
383 return(CPDIRECT_MEM_NULL);
384 }
385 cp_mem = (CPodeMem) cpode_mem;
386
387 if (lmemP == NULL) {
388 cpProcessError(cp_mem, CPDIRECT_LMEM_NULL, "CPDIRECT", "CPDlsProjGetNumFctEvals", MSGD_LMEM_NULL);
389 return(CPDIRECT_LMEM_NULL);
390 }
391 cpdlsP_mem = (CPDlsProjMem) lmemP;
392
393 *ncevalsLS = nceDQ;
394
395 return(CPDIRECT_SUCCESS);
396 }
397
398 /*
399 * =================================================================
400 * DENSE DQ JACOBIAN APPROXIMATIONS FOR IMPLICIT INTEGRATION
401 * =================================================================
402 */
403
404 /*
405 * -----------------------------------------------------------------
406 * cpDlsDenseDQJacExpl
407 * -----------------------------------------------------------------
408 * This routine generates a dense difference quotient approximation to
409 * the Jacobian of f(t,y). It assumes that a dense matrix of type
410 * DlsMat is stored column-wise, and that elements within each column
411 * are contiguous. The address of the jth column of J is obtained via
412 * the macro DENSE_COL and this pointer is associated with an N_Vector
413 * using the N_VGetArrayPointer/N_VSetArrayPointer functions.
414 * Finally, the actual computation of the jth column of the Jacobian is
415 * done with a call to N_VLinearSum.
416 * -----------------------------------------------------------------
417 */
418
cpDlsDenseDQJacExpl(int N,realtype t,N_Vector y,N_Vector fy,DlsMat Jac,void * jac_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)419 int cpDlsDenseDQJacExpl(int N, realtype t,
420 N_Vector y, N_Vector fy,
421 DlsMat Jac, void *jac_data,
422 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
423 {
424 realtype fnorm, minInc, inc, inc_inv, yjsaved, srur;
425 realtype *tmp2_data, *y_data, *ewt_data;
426 N_Vector ftemp, jthCol;
427 int j;
428 int retval = 0;
429
430 CPodeMem cp_mem;
431 CPDlsMem cpdls_mem;
432
433 /* jac_data points to cpode_mem */
434 cp_mem = (CPodeMem) jac_data;
435 cpdls_mem = (CPDlsMem) lmem;
436
437 /* Save pointer to the array in tmp2 */
438 tmp2_data = N_VGetArrayPointer(tmp2);
439
440 /* Rename work vectors for readibility */
441 ftemp = tmp1;
442 jthCol = tmp2;
443
444 /* Obtain pointers to the data for ewt, y */
445 ewt_data = N_VGetArrayPointer(ewt);
446 y_data = N_VGetArrayPointer(y);
447
448 /* Set minimum increment based on uround and norm of f */
449 srur = RSqrt(uround);
450 fnorm = N_VWrmsNorm(fy, ewt);
451 minInc = (fnorm != ZERO) ?
452 (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
453
454 for (j = 0; j < N; j++) {
455
456 /* Generate the jth col of J(tn,y) */
457
458 N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol);
459
460 yjsaved = y_data[j];
461 inc = MAX(srur*ABS(yjsaved), minInc/ewt_data[j]);
462 y_data[j] += inc;
463
464 retval = fe(t, y, ftemp, f_data);
465 nfeDQ++;
466 if (retval != 0) break;
467
468 y_data[j] = yjsaved;
469
470 inc_inv = ONE/inc;
471 N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);
472
473 DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);
474 }
475
476 /* Restore original array pointer in tmp2 */
477 N_VSetArrayPointer(tmp2_data, tmp2);
478
479 return(retval);
480 }
481
482 /*
483 * -----------------------------------------------------------------
484 * cpDlsDenseDQJacImpl
485 * -----------------------------------------------------------------
486 * This routine generates a dense difference quotient approximation to
487 * the Jacobian F_y' + gamma*F_y. It assumes that a dense matrix of type
488 * DlsMat is stored column-wise, and that elements within each column
489 * are contiguous. The address of the jth column of J is obtained via
490 * the macro DENSE_COL and this pointer is associated with an N_Vector
491 * using the N_VGetArrayPointer/N_VSetArrayPointer functions.
492 * Finally, the actual computation of the jth column of the Jacobian is
493 * done with a call to N_VLinearSum.
494 * -----------------------------------------------------------------
495 */
cpDlsDenseDQJacImpl(int N,realtype t,realtype gm,N_Vector y,N_Vector yp,N_Vector r,DlsMat Jac,void * jac_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)496 int cpDlsDenseDQJacImpl(int N, realtype t, realtype gm,
497 N_Vector y, N_Vector yp, N_Vector r,
498 DlsMat Jac, void *jac_data,
499 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
500 {
501 realtype inc, inc_inv, yj, ypj, srur;
502 realtype *tmp2_data, *y_data, *yp_data, *ewt_data;
503 N_Vector ftemp, jthCol;
504 int j;
505 int retval = 0;
506
507 CPodeMem cp_mem;
508 CPDlsMem cpdls_mem;
509
510 /* jac_data points to cpode_mem */
511 cp_mem = (CPodeMem) jac_data;
512 cpdls_mem = (CPDlsMem) lmem;
513
514 /* Save pointer to the array in tmp2 */
515 tmp2_data = N_VGetArrayPointer(tmp2);
516
517 /* Rename work vectors for readibility */
518 ftemp = tmp1;
519 jthCol = tmp2;
520
521 /* Obtain pointers to the data for ewt, y, and yp */
522 ewt_data = N_VGetArrayPointer(ewt);
523 y_data = N_VGetArrayPointer(y);
524 yp_data = N_VGetArrayPointer(yp);
525
526 /* Set minimum increment based on uround and norm of f */
527 srur = RSqrt(uround);
528
529 /* Generate each column of the Jacobian M = F_y' + gamma * F_y
530 as delta(F)/delta(y_j). */
531 for (j = 0; j < N; j++) {
532
533 /* Set data address of jthCol, and save y_j and yp_j values. */
534 N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol);
535 yj = y_data[j];
536 ypj = yp_data[j];
537
538 /* Set increment inc to y_j based on sqrt(uround)*abs(y_j), with
539 adjustments using yp_j and ewt_j if this is small, and a further
540 adjustment to give it the same sign as h*yp_j. */
541 inc = MAX( srur * MAX( ABS(yj), ABS(h*ypj) ) , ONE/ewt_data[j] );
542 if (h*ypj < ZERO) inc = -inc;
543 inc = (yj + inc) - yj;
544
545 /* Increment y_j and yp_j, call res, and break on error return. */
546 y_data[j] += gamma*inc;
547 yp_data[j] += inc;
548 retval = fi(t, y, yp, ftemp, f_data);
549 nfeDQ++;
550 if (retval != 0) break;
551
552 /* Generate the jth col of J(tn,y) */
553 inc_inv = ONE/inc;
554 N_VLinearSum(inc_inv, ftemp, -inc_inv, r, jthCol);
555
556 DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);
557
558 /* Reset y_j, yp_j */
559 y_data[j] = yj;
560 yp_data[j] = ypj;
561 }
562
563 /* Restore original array pointer in tmp2 */
564 N_VSetArrayPointer(tmp2_data, tmp2);
565
566 return(retval);
567 }
568
569 /*
570 * =================================================================
571 * BAND DQ JACOBIAN APPROXIMATIONS FOR IMPLICIT INTEGRATION
572 * =================================================================
573 */
574
575 /*
576 * -----------------------------------------------------------------
577 * cpDlsBandDQJacExpl
578 * -----------------------------------------------------------------
579 * This routine generates a banded difference quotient approximation to
580 * the Jacobian of f(t,y). It assumes that a band matrix of type
581 * DlsMat is stored column-wise, and that elements within each column
582 * are contiguous. This makes it possible to get the address of a column
583 * of J via the macro BAND_COL and to write a simple for loop to set
584 * each of the elements of a column in succession.
585 * -----------------------------------------------------------------
586 */
587
cpDlsBandDQJacExpl(int N,int mupper,int mlower,realtype t,N_Vector y,N_Vector fy,DlsMat Jac,void * jac_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)588 int cpDlsBandDQJacExpl(int N, int mupper, int mlower,
589 realtype t, N_Vector y, N_Vector fy,
590 DlsMat Jac, void *jac_data,
591 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
592 {
593 N_Vector ftemp, ytemp;
594 realtype fnorm, minInc, inc, inc_inv, srur;
595 realtype *col_j, *ewt_data, *fy_data, *ftemp_data, *y_data, *ytemp_data;
596 int group, i, j, width, ngroups, i1, i2;
597 int retval = 0;
598
599 CPodeMem cp_mem;
600 CPDlsMem cpdls_mem;
601
602 /* jac_dat points to cpode_mem */
603 cp_mem = (CPodeMem) jac_data;
604 cpdls_mem = (CPDlsMem) lmem;
605
606 /* Rename work vectors for use as temporary values of y and f */
607 ftemp = tmp1;
608 ytemp = tmp2;
609
610 /* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp */
611 ewt_data = N_VGetArrayPointer(ewt);
612 fy_data = N_VGetArrayPointer(fy);
613 ftemp_data = N_VGetArrayPointer(ftemp);
614 y_data = N_VGetArrayPointer(y);
615 ytemp_data = N_VGetArrayPointer(ytemp);
616
617 /* Load ytemp with y = predicted y vector */
618 N_VScale(ONE, y, ytemp);
619
620 /* Set minimum increment based on uround and norm of f */
621 srur = RSqrt(uround);
622 fnorm = N_VWrmsNorm(fy, ewt);
623 minInc = (fnorm != ZERO) ?
624 (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
625
626 /* Set bandwidth and number of column groups for band differencing */
627 width = mlower + mupper + 1;
628 ngroups = MIN(width, N);
629
630 /* Loop over column groups. */
631 for (group=1; group <= ngroups; group++) {
632
633 /* Increment all y_j in group */
634 for(j=group-1; j < N; j+=width) {
635 inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
636 ytemp_data[j] += inc;
637 }
638
639 /* Evaluate f with incremented y */
640
641 retval = fe(tn, ytemp, ftemp, f_data);
642 nfeDQ++;
643 if (retval != 0) break;
644
645 /* Restore ytemp, then form and load difference quotients */
646 for (j=group-1; j < N; j+=width) {
647 ytemp_data[j] = y_data[j];
648 col_j = BAND_COL(Jac,j);
649 inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
650 inc_inv = ONE/inc;
651 i1 = MAX(0, j-mupper);
652 i2 = MIN(j+mlower, N-1);
653 for (i=i1; i <= i2; i++)
654 BAND_COL_ELEM(col_j,i,j) = inc_inv * (ftemp_data[i] - fy_data[i]);
655 }
656 }
657
658 return(retval);
659 }
660
661 /*
662 * -----------------------------------------------------------------
663 * cpDlsBandDQJacImpl
664 * -----------------------------------------------------------------
665 */
666
cpDlsBandDQJacImpl(int N,int mupper,int mlower,realtype t,realtype gm,N_Vector y,N_Vector yp,N_Vector r,DlsMat Jac,void * jac_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)667 int cpDlsBandDQJacImpl(int N, int mupper, int mlower,
668 realtype t, realtype gm,
669 N_Vector y, N_Vector yp, N_Vector r,
670 DlsMat Jac, void *jac_data,
671 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
672 {
673 N_Vector ftemp, ytemp, yptemp;
674 realtype inc, inc_inv, yj, ypj, srur, ewtj;
675 realtype *y_data, *yp_data, *ewt_data;
676 realtype *ytemp_data, *yptemp_data, *ftemp_data, *r_data, *col_j;
677 int i, j, i1, i2, width, group, ngroups;
678 int retval = 0;
679
680 CPodeMem cp_mem;
681 CPDlsMem cpdls_mem;
682
683 /* jac_data points to cpode_mem */
684 cp_mem = (CPodeMem) jac_data;
685 cpdls_mem = (CPDlsMem) lmem;
686
687 ftemp = tmp1; /* Rename work vector for use as the perturbed residual. */
688 ytemp = tmp2; /* Rename work vector for use as a temporary for yy. */
689 yptemp= tmp3; /* Rename work vector for use as a temporary for yp. */
690
691 /* Obtain pointers to the data for all vectors used. */
692 ewt_data = N_VGetArrayPointer(ewt);
693 r_data = N_VGetArrayPointer(r);
694 y_data = N_VGetArrayPointer(y);
695 yp_data = N_VGetArrayPointer(yp);
696 ftemp_data = N_VGetArrayPointer(ftemp);
697 ytemp_data = N_VGetArrayPointer(ytemp);
698 yptemp_data = N_VGetArrayPointer(yptemp);
699
700 /* Initialize ytemp and yptemp. */
701 N_VScale(ONE, y, ytemp);
702 N_VScale(ONE, yp, yptemp);
703
704 /* Compute miscellaneous values for the Jacobian computation. */
705 srur = RSqrt(uround);
706 width = mlower + mupper + 1;
707 ngroups = MIN(width, N);
708
709 /* Loop over column groups. */
710 for (group=1; group <= ngroups; group++) {
711
712 /* Increment all y[j] and yp[j] for j in this group. */
713 for (j=group-1; j<N; j+=width) {
714 yj = y_data[j];
715 ypj = yp_data[j];
716 ewtj = ewt_data[j];
717
718 /* Set increment inc to yj based on sqrt(uround)*abs(yj), with
719 adjustments using ypj and ewtj if this is small, and a further
720 adjustment to give it the same sign as h*ypj. */
721 inc = MAX( srur * MAX( ABS(yj), ABS(h*ypj) ) , ONE/ewtj );
722
723 if (h*ypj < ZERO) inc = -inc;
724 inc = (yj + inc) - yj;
725
726 /* Increment yj and ypj. */
727 ytemp_data[j] += gamma*inc;
728 yptemp_data[j] += inc;
729 }
730
731 /* Call ODE fct. with incremented arguments. */
732 retval = fi(tn, ytemp, yptemp, ftemp, f_data);
733 nfeDQ++;
734 if (retval != 0) break;
735
736 /* Loop over the indices j in this group again. */
737 for (j=group-1; j<N; j+=width) {
738
739 /* Reset ytemp and yptemp components that were perturbed. */
740 yj = ytemp_data[j] = y_data[j];
741 ypj = yptemp_data[j] = yp_data[j];
742 col_j = BAND_COL(Jac, j);
743 ewtj = ewt_data[j];
744
745 /* Set increment inc exactly as above. */
746 inc = MAX( srur * MAX( ABS(yj), ABS(h*ypj) ) , ONE/ewtj );
747 if (h*ypj < ZERO) inc = -inc;
748 inc = (yj + inc) - yj;
749
750 /* Load the difference quotient Jacobian elements for column j. */
751 inc_inv = ONE/inc;
752 i1 = MAX(0, j-mupper);
753 i2 = MIN(j+mlower,N-1);
754 for (i=i1; i<=i2; i++)
755 BAND_COL_ELEM(col_j,i,j) = inc_inv*(ftemp_data[i]-r_data[i]);
756 }
757
758 }
759
760 return(retval);
761 }
762
763 /*
764 * =================================================================
765 * DENSE DQ JACOBIAN APPROXIMATIONS FOR PROJECTION
766 * =================================================================
767 */
768
769 /*
770 * -----------------------------------------------------------------
771 * cpDlsDenseProjDQJac
772 * -----------------------------------------------------------------
773 * This routine generates a dense difference quotient approximation
774 * to the transpose of the Jacobian of c(t,y). It loads it into a
775 * dense matrix of type DlsMat stored column-wise with elements
776 * within each column contiguous. The address of the jth column of
777 * J is obtained via the macro DENSE_COL and this pointer is
778 * associated with an N_Vector using the N_VGetArrayPointer and
779 * N_VSetArrayPointer functions.
780 * Finally, the actual computation of the jth column of the Jacobian
781 * transposed is done with a call to N_VLinearSum.
782 * -----------------------------------------------------------------
783 */
cpDlsDenseProjDQJac(int Nc,int Ny,realtype t,N_Vector y,N_Vector cy,DlsMat Jac,void * jac_data,N_Vector c_tmp1,N_Vector c_tmp2)784 int cpDlsDenseProjDQJac(int Nc, int Ny, realtype t,
785 N_Vector y, N_Vector cy,
786 DlsMat Jac, void *jac_data,
787 N_Vector c_tmp1, N_Vector c_tmp2)
788 {
789 realtype inc, inc_inv, yj, srur;
790 realtype *y_data, *ewt_data, *jthCol_data;
791 N_Vector ctemp, jthCol;
792 int i, j;
793 int retval = 0;
794
795 CPodeMem cp_mem;
796 CPDlsProjMem cpdlsP_mem;
797
798 /* jac_data points to cpode_mem */
799 cp_mem = (CPodeMem) jac_data;
800 cpdlsP_mem = (CPDlsProjMem) lmemP;
801
802 /* Rename work vectors for readibility */
803 ctemp = c_tmp1;
804 jthCol = c_tmp2;
805
806 /* Obtain pointers to the data for ewt and y */
807 ewt_data = N_VGetArrayPointer(ewt);
808 y_data = N_VGetArrayPointer(y);
809
810 /* Obtain pointer to the data for jthCol */
811 jthCol_data = N_VGetArrayPointer(jthCol);
812
813 /* Set minimum increment based on uround and norm of f */
814 srur = RSqrt(uround);
815
816 /* Generate each column of the Jacobian G = dc/dy as delta(c)/delta(y_j). */
817 for (j = 0; j < Ny; j++) {
818
819 /* Save the y_j values. */
820 yj = y_data[j];
821
822 /* Set increment inc to y_j based on sqrt(uround)*abs(y_j),
823 with an adjustment using ewt_j if this is small */
824 inc = MAX( srur * ABS(yj) , ONE/ewt_data[j] );
825 inc = (yj + inc) - yj;
826
827 /* Increment y_j, call cfun, and break on error return. */
828 y_data[j] += inc;
829 retval = cfun(t, y, ctemp, c_data);
830 nceDQ++;
831 if (retval != 0) break;
832
833 /* Generate the jth col of G(tn,y) */
834 inc_inv = ONE/inc;
835 N_VLinearSum(inc_inv, ctemp, -inc_inv, cy, jthCol);
836
837 /* Copy the j-th column of G into the j-th row of Jac */
838 for (i = 0; i < Nc ; i++) {
839 DENSE_ELEM(Jac,j,i) = jthCol_data[i];
840 }
841
842 /* Reset y_j */
843 y_data[j] = yj;
844 }
845
846 return(retval);
847
848 }
849
850
851