1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 1.6 $
4 * $Date: 2010/12/01 22:35:26 $
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 an IDADLS linear solver.
14 * -----------------------------------------------------------------
15 */
16
17 /*
18 * =================================================================
19 * IMPORTED HEADER FILES
20 * =================================================================
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "ida_impl.h"
27 #include "ida_direct_impl.h"
28 #include <sundials/sundials_math.h>
29
30 /*
31 * =================================================================
32 * FUNCTION SPECIFIC CONSTANTS
33 * =================================================================
34 */
35
36 #define ZERO RCONST(0.0)
37 #define ONE RCONST(1.0)
38 #define TWO RCONST(2.0)
39
40 /*
41 * =================================================================
42 * READIBILITY REPLACEMENTS
43 * =================================================================
44 */
45
46 #define res (IDA_mem->ida_res)
47 #define user_data (IDA_mem->ida_user_data)
48 #define uround (IDA_mem->ida_uround)
49 #define nst (IDA_mem->ida_nst)
50 #define tn (IDA_mem->ida_tn)
51 #define hh (IDA_mem->ida_hh)
52 #define cj (IDA_mem->ida_cj)
53 #define cjratio (IDA_mem->ida_cjratio)
54 #define ewt (IDA_mem->ida_ewt)
55 #define constraints (IDA_mem->ida_constraints)
56
57 #define linit (IDA_mem->ida_linit)
58 #define lsetup (IDA_mem->ida_lsetup)
59 #define lsolve (IDA_mem->ida_lsolve)
60 #define lfree (IDA_mem->ida_lfree)
61 #define lperf (IDA_mem->ida_lperf)
62 #define lmem (IDA_mem->ida_lmem)
63 #define tempv (IDA_mem->ida_tempv1)
64 #define setupNonNull (IDA_mem->ida_setupNonNull)
65
66 #define mtype (idadls_mem->d_type)
67 #define n (idadls_mem->d_n)
68 #define ml (idadls_mem->d_ml)
69 #define mu (idadls_mem->d_mu)
70 #define smu (idadls_mem->d_smu)
71 #define jacDQ (idadls_mem->d_jacDQ)
72 #define djac (idadls_mem->d_djac)
73 #define bjac (idadls_mem->d_bjac)
74 #define M (idadls_mem->d_J)
75 #define pivots (idadls_mem->d_pivots)
76 #define nje (idadls_mem->d_nje)
77 #define nreDQ (idadls_mem->d_nreDQ)
78 #define last_flag (idadls_mem->d_last_flag)
79
80 /*
81 * =================================================================
82 * EXPORTED FUNCTIONS FOR IMPLICIT INTEGRATION
83 * =================================================================
84 */
85
86 /*
87 * IDADlsSetDenseJacFn specifies the dense Jacobian function.
88 */
IDADlsSetDenseJacFn(void * ida_mem,IDADlsDenseJacFn jac)89 int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn jac)
90 {
91 IDAMem IDA_mem;
92 IDADlsMem idadls_mem;
93
94 /* Return immediately if ida_mem is NULL */
95 if (ida_mem == NULL)
96 {
97 IDAProcessError(NULL, IDADLS_MEM_NULL, "IDADLS", "IDADlsSetDenseJacFn", MSGD_IDAMEM_NULL);
98 return(IDADLS_MEM_NULL);
99 }
100 IDA_mem = (IDAMem) ida_mem;
101
102 if (lmem == NULL)
103 {
104 IDAProcessError(IDA_mem, IDADLS_LMEM_NULL, "IDADLS", "IDADlsSetDenseJacFn", MSGD_LMEM_NULL);
105 return(IDADLS_LMEM_NULL);
106 }
107 idadls_mem = (IDADlsMem) lmem;
108
109 if (jac != NULL)
110 {
111 jacDQ = FALSE;
112 djac = jac;
113 }
114 else
115 {
116 jacDQ = TRUE;
117 }
118
119 return(IDADLS_SUCCESS);
120 }
121
122 /*
123 * IDADlsSetBandJacFn specifies the band Jacobian function.
124 */
IDADlsSetBandJacFn(void * ida_mem,IDADlsBandJacFn jac)125 int IDADlsSetBandJacFn(void *ida_mem, IDADlsBandJacFn jac)
126 {
127 IDAMem IDA_mem;
128 IDADlsMem idadls_mem;
129
130 /* Return immediately if ida_mem is NULL */
131 if (ida_mem == NULL)
132 {
133 IDAProcessError(NULL, IDADLS_MEM_NULL, "IDADLS", "IDADlsSetBandJacFn", MSGD_IDAMEM_NULL);
134 return(IDADLS_MEM_NULL);
135 }
136 IDA_mem = (IDAMem) ida_mem;
137
138 if (lmem == NULL)
139 {
140 IDAProcessError(IDA_mem, IDADLS_LMEM_NULL, "IDADLS", "IDADlsSetBandJacFn", MSGD_LMEM_NULL);
141 return(IDADLS_LMEM_NULL);
142 }
143 idadls_mem = (IDADlsMem) lmem;
144
145 if (jac != NULL)
146 {
147 jacDQ = FALSE;
148 bjac = jac;
149 }
150 else
151 {
152 jacDQ = TRUE;
153 }
154
155 return(IDADLS_SUCCESS);
156 }
157
158 /*
159 * IDADlsGetWorkSpace returns the length of workspace allocated for the
160 * IDALAPACK linear solver.
161 */
IDADlsGetWorkSpace(void * ida_mem,long int * lenrwLS,long int * leniwLS)162 int IDADlsGetWorkSpace(void *ida_mem, long int *lenrwLS, long int *leniwLS)
163 {
164 IDAMem IDA_mem;
165 IDADlsMem idadls_mem;
166
167 /* Return immediately if ida_mem is NULL */
168 if (ida_mem == NULL)
169 {
170 IDAProcessError(NULL, IDADLS_MEM_NULL, "IDADLS", "IDADlsGetWorkSpace", MSGD_IDAMEM_NULL);
171 return(IDADLS_MEM_NULL);
172 }
173 IDA_mem = (IDAMem) ida_mem;
174
175 if (lmem == NULL)
176 {
177 IDAProcessError(IDA_mem, IDADLS_LMEM_NULL, "IDADLS", "IDADlsGetWorkSpace", MSGD_LMEM_NULL);
178 return(IDADLS_LMEM_NULL);
179 }
180 idadls_mem = (IDADlsMem) lmem;
181
182 if (mtype == SUNDIALS_DENSE)
183 {
184 *lenrwLS = n * n;
185 *leniwLS = n;
186 }
187 else if (mtype == SUNDIALS_BAND)
188 {
189 *lenrwLS = n * (smu + ml + 1);
190 *leniwLS = n;
191 }
192
193 return(IDADLS_SUCCESS);
194 }
195
196 /*
197 * IDADlsGetNumJacEvals returns the number of Jacobian evaluations.
198 */
IDADlsGetNumJacEvals(void * ida_mem,long int * njevals)199 int IDADlsGetNumJacEvals(void *ida_mem, long int *njevals)
200 {
201 IDAMem IDA_mem;
202 IDADlsMem idadls_mem;
203
204 /* Return immediately if ida_mem is NULL */
205 if (ida_mem == NULL)
206 {
207 IDAProcessError(NULL, IDADLS_MEM_NULL, "IDADLS", "IDADlsGetNumJacEvals", MSGD_IDAMEM_NULL);
208 return(IDADLS_MEM_NULL);
209 }
210 IDA_mem = (IDAMem) ida_mem;
211
212 if (lmem == NULL)
213 {
214 IDAProcessError(IDA_mem, IDADLS_LMEM_NULL, "IDADLS", "IDADlsGetNumJacEvals", MSGD_LMEM_NULL);
215 return(IDADLS_LMEM_NULL);
216 }
217 idadls_mem = (IDADlsMem) lmem;
218
219 *njevals = nje;
220
221 return(IDADLS_SUCCESS);
222 }
223
224 /*
225 * IDADlsGetNumResEvals returns the number of calls to the DAE function
226 * needed for the DQ Jacobian approximation.
227 */
IDADlsGetNumResEvals(void * ida_mem,long int * nrevalsLS)228 int IDADlsGetNumResEvals(void *ida_mem, long int *nrevalsLS)
229 {
230 IDAMem IDA_mem;
231 IDADlsMem idadls_mem;
232
233 /* Return immediately if ida_mem is NULL */
234 if (ida_mem == NULL)
235 {
236 IDAProcessError(NULL, IDADLS_MEM_NULL, "IDADLS", "IDADlsGetNumFctEvals", MSGD_IDAMEM_NULL);
237 return(IDADLS_MEM_NULL);
238 }
239 IDA_mem = (IDAMem) ida_mem;
240
241 if (lmem == NULL)
242 {
243 IDAProcessError(IDA_mem, IDADLS_LMEM_NULL, "IDADLS", "IDADlsGetNumFctEvals", MSGD_LMEM_NULL);
244 return(IDADLS_LMEM_NULL);
245 }
246 idadls_mem = (IDADlsMem) lmem;
247
248 *nrevalsLS = nreDQ;
249
250 return(IDADLS_SUCCESS);
251 }
252
253 /*
254 * IDADlsGetReturnFlagName returns the name associated with a IDALAPACK
255 * return value.
256 */
IDADlsGetReturnFlagName(long int flag)257 char *IDADlsGetReturnFlagName(long int flag)
258 {
259 char *name;
260
261 name = (char *)malloc(30 * sizeof(char));
262
263 switch (flag)
264 {
265 case IDADLS_SUCCESS:
266 sprintf(name, "IDADLS_SUCCESS");
267 break;
268 case IDADLS_MEM_NULL:
269 sprintf(name, "IDADLS_MEM_NULL");
270 break;
271 case IDADLS_LMEM_NULL:
272 sprintf(name, "IDADLS_LMEM_NULL");
273 break;
274 case IDADLS_ILL_INPUT:
275 sprintf(name, "IDADLS_ILL_INPUT");
276 break;
277 case IDADLS_MEM_FAIL:
278 sprintf(name, "IDADLS_MEM_FAIL");
279 break;
280 case IDADLS_JACFUNC_UNRECVR:
281 sprintf(name, "IDADLS_JACFUNC_UNRECVR");
282 break;
283 case IDADLS_JACFUNC_RECVR:
284 sprintf(name, "IDADLS_JACFUNC_RECVR");
285 break;
286 default:
287 sprintf(name, "NONE");
288 }
289
290 return(name);
291 }
292
293 /*
294 * IDADlsGetLastFlag returns the last flag set in a IDALAPACK function.
295 */
IDADlsGetLastFlag(void * ida_mem,long int * flag)296 int IDADlsGetLastFlag(void *ida_mem, long int *flag)
297 {
298 IDAMem IDA_mem;
299 IDADlsMem idadls_mem;
300
301 /* Return immediately if ida_mem is NULL */
302 if (ida_mem == NULL)
303 {
304 IDAProcessError(NULL, IDADLS_MEM_NULL, "IDADLS", "IDADlsGetLastFlag", MSGD_IDAMEM_NULL);
305 return(IDADLS_MEM_NULL);
306 }
307 IDA_mem = (IDAMem) ida_mem;
308
309 if (lmem == NULL)
310 {
311 IDAProcessError(IDA_mem, IDADLS_LMEM_NULL, "IDADLS", "IDADlsGetLastFlag", MSGD_LMEM_NULL);
312 return(IDADLS_LMEM_NULL);
313 }
314 idadls_mem = (IDADlsMem) lmem;
315
316 *flag = last_flag;
317
318 return(IDADLS_SUCCESS);
319 }
320
321 /*
322 * =================================================================
323 * DQ JACOBIAN APPROXIMATIONS
324 * =================================================================
325 */
326
327 /*
328 * -----------------------------------------------------------------
329 * idaDlsDenseDQJac
330 * -----------------------------------------------------------------
331 * This routine generates a dense difference quotient approximation to
332 * the Jacobian F_y + c_j*F_y'. It assumes that a dense matrix of type
333 * DlsMat is stored column-wise, and that elements within each column
334 * are contiguous. The address of the jth column of J is obtained via
335 * the macro LAPACK_DENSE_COL and this pointer is associated with an N_Vector
336 * using the N_VGetArrayPointer/N_VSetArrayPointer functions.
337 * Finally, the actual computation of the jth column of the Jacobian is
338 * done with a call to N_VLinearSum.
339 * -----------------------------------------------------------------
340 */
idaDlsDenseDQJac(long int N,realtype tt,realtype c_j,N_Vector yy,N_Vector yp,N_Vector rr,DlsMat Jac,void * data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)341 int idaDlsDenseDQJac(long int N, realtype tt, realtype c_j,
342 N_Vector yy, N_Vector yp, N_Vector rr,
343 DlsMat Jac, void *data,
344 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
345 {
346 realtype inc, inc_inv, yj, ypj, srur, conj;
347 realtype *tmp2_data, *y_data, *yp_data, *ewt_data, *cns_data = NULL;
348 N_Vector rtemp, jthCol;
349 long int j;
350 int retval = 0;
351
352 IDAMem IDA_mem;
353 IDADlsMem idadls_mem;
354
355 /* data points to IDA_mem */
356 IDA_mem = (IDAMem) data;
357 idadls_mem = (IDADlsMem) lmem;
358
359 /* Save pointer to the array in tmp2 */
360 tmp2_data = N_VGetArrayPointer(tmp2);
361
362 /* Rename work vectors for readibility */
363 rtemp = tmp1;
364 jthCol = tmp2;
365
366 /* Obtain pointers to the data for ewt, yy, yp. */
367 ewt_data = N_VGetArrayPointer(ewt);
368 y_data = N_VGetArrayPointer(yy);
369 yp_data = N_VGetArrayPointer(yp);
370 if (constraints != NULL)
371 {
372 cns_data = N_VGetArrayPointer(constraints);
373 }
374
375 srur = RSqrt(uround);
376
377 for (j = 0; j < N; j++)
378 {
379
380 /* Generate the jth col of J(tt,yy,yp) as delta(F)/delta(y_j). */
381
382 /* Set data address of jthCol, and save y_j and yp_j values. */
383 N_VSetArrayPointer(DENSE_COL(Jac, j), jthCol);
384 yj = y_data[j];
385 ypj = yp_data[j];
386
387 /* Set increment inc to y_j based on sqrt(uround)*abs(y_j), with
388 adjustments using yp_j and ewt_j if this is small, and a further
389 adjustment to give it the same sign as hh*yp_j. */
390
391 inc = MAX( srur * MAX( ABS(yj), ABS(hh * ypj) ) , ONE / ewt_data[j] );
392
393 if (hh * ypj < ZERO)
394 {
395 inc = -inc;
396 }
397 inc = (yj + inc) - yj;
398
399 /* Adjust sign(inc) again if y_j has an inequality constraint. */
400 if (constraints != NULL)
401 {
402 conj = cns_data[j];
403 if (ABS(conj) == ONE)
404 {
405 if ((yj + inc)*conj < ZERO)
406 {
407 inc = -inc;
408 }
409 }
410 else if (ABS(conj) == TWO)
411 {
412 if ((yj + inc)*conj <= ZERO)
413 {
414 inc = -inc;
415 }
416 }
417 }
418
419 /* Increment y_j and yp_j, call res, and break on error return. */
420 y_data[j] += inc;
421 yp_data[j] += c_j * inc;
422
423 retval = res(tt, yy, yp, rtemp, user_data);
424 nreDQ++;
425 if (retval != 0)
426 {
427 break;
428 }
429
430 /* Construct difference quotient in jthCol */
431 inc_inv = ONE / inc;
432 N_VLinearSum(inc_inv, rtemp, -inc_inv, rr, jthCol);
433
434 DENSE_COL(Jac, j) = N_VGetArrayPointer(jthCol);
435
436 /* reset y_j, yp_j */
437 y_data[j] = yj;
438 yp_data[j] = ypj;
439 }
440
441 /* Restore original array pointer in tmp2 */
442 N_VSetArrayPointer(tmp2_data, tmp2);
443
444 return(retval);
445
446 }
447
448 /*
449 * -----------------------------------------------------------------
450 * idaDlsBandDQJac
451 * -----------------------------------------------------------------
452 * This routine generates a banded difference quotient approximation JJ
453 * to the DAE system Jacobian J. It assumes that a band matrix of type
454 * BandMat is stored column-wise, and that elements within each column
455 * are contiguous. The address of the jth column of JJ is obtained via
456 * the macros BAND_COL and BAND_COL_ELEM. The columns of the Jacobian are
457 * constructed using mupper + mlower + 1 calls to the res routine, and
458 * appropriate differencing.
459 * The return value is either IDABAND_SUCCESS = 0, or the nonzero value returned
460 * by the res routine, if any.
461 */
462
idaDlsBandDQJac(long int N,long int mupper,long int mlower,realtype tt,realtype c_j,N_Vector yy,N_Vector yp,N_Vector rr,DlsMat Jac,void * data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)463 int idaDlsBandDQJac(long int N, long int mupper, long int mlower,
464 realtype tt, realtype c_j,
465 N_Vector yy, N_Vector yp, N_Vector rr,
466 DlsMat Jac, void *data,
467 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
468 {
469 realtype inc, inc_inv, yj, ypj, srur, conj, ewtj;
470 realtype *y_data, *yp_data, *ewt_data, *cns_data = NULL;
471 realtype *ytemp_data, *yptemp_data, *rtemp_data, *r_data, *col_j;
472 N_Vector rtemp, ytemp, yptemp;
473 long int i, j, i1, i2, width, ngroups, group;
474 int retval = 0;
475
476 IDAMem IDA_mem;
477 IDADlsMem idadls_mem;
478
479 /* data points to IDA_mem */
480 IDA_mem = (IDAMem) data;
481 idadls_mem = (IDADlsMem) lmem;
482
483 rtemp = tmp1; /* Rename work vector for use as the perturbed residual. */
484
485 ytemp = tmp2; /* Rename work vector for use as a temporary for yy. */
486
487
488 yptemp = tmp3; /* Rename work vector for use as a temporary for yp. */
489
490 /* Obtain pointers to the data for all eight vectors used. */
491
492 ewt_data = N_VGetArrayPointer(ewt);
493 r_data = N_VGetArrayPointer(rr);
494 y_data = N_VGetArrayPointer(yy);
495 yp_data = N_VGetArrayPointer(yp);
496
497 rtemp_data = N_VGetArrayPointer(rtemp);
498 ytemp_data = N_VGetArrayPointer(ytemp);
499 yptemp_data = N_VGetArrayPointer(yptemp);
500
501 if (constraints != NULL)
502 {
503 cns_data = N_VGetArrayPointer(constraints);
504 }
505
506 /* Initialize ytemp and yptemp. */
507
508 N_VScale(ONE, yy, ytemp);
509 N_VScale(ONE, yp, yptemp);
510
511 /* Compute miscellaneous values for the Jacobian computation. */
512
513 srur = RSqrt(uround);
514 width = mlower + mupper + 1;
515 ngroups = MIN(width, N);
516
517 /* Loop over column groups. */
518 for (group = 1; group <= ngroups; group++)
519 {
520
521 /* Increment all yy[j] and yp[j] for j in this group. */
522
523 for (j = group - 1; j < N; j += width)
524 {
525 yj = y_data[j];
526 ypj = yp_data[j];
527 ewtj = ewt_data[j];
528
529 /* Set increment inc to yj based on sqrt(uround)*abs(yj), with
530 adjustments using ypj and ewtj if this is small, and a further
531 adjustment to give it the same sign as hh*ypj. */
532
533 inc = MAX( srur * MAX( ABS(yj), ABS(hh * ypj) ) , ONE / ewtj );
534
535 if (hh * ypj < ZERO)
536 {
537 inc = -inc;
538 }
539 inc = (yj + inc) - yj;
540
541 /* Adjust sign(inc) again if yj has an inequality constraint. */
542
543 if (constraints != NULL)
544 {
545 conj = cns_data[j];
546 if (ABS(conj) == ONE)
547 {
548 if ((yj + inc)*conj < ZERO)
549 {
550 inc = -inc;
551 }
552 }
553 else if (ABS(conj) == TWO)
554 {
555 if ((yj + inc)*conj <= ZERO)
556 {
557 inc = -inc;
558 }
559 }
560 }
561
562 /* Increment yj and ypj. */
563
564 ytemp_data[j] += inc;
565 yptemp_data[j] += cj * inc;
566 }
567
568 /* Call res routine with incremented arguments. */
569
570 retval = res(tt, ytemp, yptemp, rtemp, user_data);
571 nreDQ++;
572 if (retval != 0)
573 {
574 break;
575 }
576
577 /* Loop over the indices j in this group again. */
578
579 for (j = group - 1; j < N; j += width)
580 {
581
582 /* Reset ytemp and yptemp components that were perturbed. */
583
584 yj = ytemp_data[j] = y_data[j];
585 ypj = yptemp_data[j] = yp_data[j];
586 col_j = BAND_COL(Jac, j);
587 ewtj = ewt_data[j];
588
589 /* Set increment inc exactly as above. */
590
591 inc = MAX( srur * MAX( ABS(yj), ABS(hh * ypj) ) , ONE / ewtj );
592 if (hh * ypj < ZERO)
593 {
594 inc = -inc;
595 }
596 inc = (yj + inc) - yj;
597 if (constraints != NULL)
598 {
599 conj = cns_data[j];
600 if (ABS(conj) == ONE)
601 {
602 if ((yj + inc)*conj < ZERO)
603 {
604 inc = -inc;
605 }
606 }
607 else if (ABS(conj) == TWO)
608 {
609 if ((yj + inc)*conj <= ZERO)
610 {
611 inc = -inc;
612 }
613 }
614 }
615
616 /* Load the difference quotient Jacobian elements for column j. */
617
618 inc_inv = ONE / inc;
619 i1 = MAX(0, j - mupper);
620 i2 = MIN(j + mlower, N - 1);
621
622 for (i = i1; i <= i2; i++)
623 {
624 BAND_COL_ELEM(col_j, i, j) = inc_inv * (rtemp_data[i] - r_data[i]);
625 }
626 }
627
628 }
629
630 return(retval);
631
632 }
633