1 /*-----------------------------------------------------------------
2 * Programmer(s): Daniel R. Reynolds @ SMU
3 * Aaron Collier and Radu Serban @ LLNL
4 *-----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2021, Lawrence Livermore National Security
7 * and Southern Methodist University.
8 * All rights reserved.
9 *
10 * See the top-level LICENSE and NOTICE files for details.
11 *
12 * SPDX-License-Identifier: BSD-3-Clause
13 * SUNDIALS Copyright End
14 *-----------------------------------------------------------------
15 * The C functions FIDAJTSetup and FIDAJtimes are to interface
16 * between the IDALS module and the user-supplied
17 * Jacobian-vector product routines FIDAJTSETUP and FIDAJTIMES.
18 * Note the use of the generic names FIDA_JTSETUP and FIDA_JTIMES
19 * below.
20 *-----------------------------------------------------------------*/
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include "fida.h" /* actual fn. names, prototypes and global vars.*/
26 #include "ida_impl.h" /* definition of IDAMem type */
27
28 #include <ida/ida_ls.h>
29
30 /*************************************************/
31
32 #ifdef __cplusplus /* wrapper to enable C++ usage */
33 extern "C" {
34 #endif
35
36 extern void FIDA_JTSETUP(realtype *T, realtype *Y, realtype *YP,
37 realtype *R, realtype *CJ, realtype *EWT,
38 realtype *H, long int *IPAR,
39 realtype *RPAR, int *IER);
40
41 extern void FIDA_JTIMES(realtype *T, realtype *Y, realtype *YP,
42 realtype *R, realtype *V, realtype *FJV,
43 realtype *CJ, realtype *EWT, realtype *H,
44 long int *IPAR, realtype *RPAR,
45 realtype *WK1, realtype *WK2, int *IER);
46
47 #ifdef __cplusplus
48 }
49 #endif
50
51 /*************************************************/
52
53 /*** DEPRECATED ***/
FIDA_SPILSSETJAC(int * flag,int * ier)54 void FIDA_SPILSSETJAC(int *flag, int *ier)
55 { FIDA_LSSETJAC(flag, ier); }
56
57
58 /* Fortran interface to C routine IDASetJacTimes; see
59 fida.h for further information */
FIDA_LSSETJAC(int * flag,int * ier)60 void FIDA_LSSETJAC(int *flag, int *ier)
61 {
62 if (*flag == 0) {
63 *ier = IDASetJacTimes(IDA_idamem, NULL, NULL);
64 } else {
65 if (F2C_IDA_ewtvec == NULL) {
66 F2C_IDA_ewtvec = N_VClone(F2C_IDA_vec);
67 if (F2C_IDA_ewtvec == NULL) {
68 *ier = -1;
69 return;
70 }
71 }
72 *ier = IDASetJacTimes(IDA_idamem, FIDAJTSetup, FIDAJtimes);
73 }
74 return;
75 }
76
77 /*************************************************/
78
79 /* C interface to user-supplied Fortran routine FIDAJTSETUP; see
80 fida.h for further information */
FIDAJTSetup(realtype t,N_Vector y,N_Vector yp,N_Vector r,realtype cj,void * user_data)81 int FIDAJTSetup(realtype t, N_Vector y, N_Vector yp,
82 N_Vector r, realtype cj, void *user_data)
83 {
84 realtype *ydata, *ypdata, *rdata, *ewtdata;
85 realtype h;
86 FIDAUserData IDA_userdata;
87 int ier = 0;
88
89 /* Initialize all pointers to NULL */
90 ydata = ypdata = rdata = ewtdata = NULL;
91
92 /* NOTE: The user-supplied routine should set ier to an
93 appropriate value, but we preset the value to zero
94 (meaning SUCCESS) so the user need only reset the
95 value if an error occurred */
96 ier = 0;
97
98 IDAGetErrWeights(IDA_idamem, F2C_IDA_ewtvec);
99 IDAGetLastStep(IDA_idamem, &h);
100 ydata = N_VGetArrayPointer(y);
101 ypdata = N_VGetArrayPointer(yp);
102 rdata = N_VGetArrayPointer(r);
103 ewtdata = N_VGetArrayPointer(F2C_IDA_ewtvec);
104
105 IDA_userdata = (FIDAUserData) user_data;
106
107 /* Call user-supplied routine */
108 FIDA_JTSETUP(&t, ydata, ypdata, rdata, &cj, ewtdata, &h,
109 IDA_userdata->ipar, IDA_userdata->rpar, &ier);
110 return(ier);
111 }
112
113 /* C interface to user-supplied Fortran routine FIDAJTIMES; see
114 fida.h for further information */
FIDAJtimes(realtype t,N_Vector yy,N_Vector yp,N_Vector rr,N_Vector v,N_Vector Jv,realtype c_j,void * user_data,N_Vector vtemp1,N_Vector vtemp2)115 int FIDAJtimes(realtype t, N_Vector yy, N_Vector yp, N_Vector rr,
116 N_Vector v, N_Vector Jv, realtype c_j,
117 void *user_data, N_Vector vtemp1, N_Vector vtemp2)
118 {
119 realtype *yy_data, *yp_data, *rr_data, *vdata, *Jvdata, *ewtdata;
120 realtype *v1data, *v2data;
121 realtype h;
122 FIDAUserData IDA_userdata;
123 int ier;
124
125 /* Initialize all pointers to NULL */
126 yy_data = yp_data = rr_data = vdata = Jvdata = ewtdata = NULL;
127
128 /* NOTE: The user-supplied routine should set ier to an
129 appropriate value, but we preset the value to zero
130 (meaning SUCCESS) so the user need only reset the
131 value if an error occurred */
132 ier = 0;
133
134 IDAGetErrWeights(IDA_idamem, F2C_IDA_ewtvec);
135 IDAGetLastStep(IDA_idamem, &h);
136
137 /* Get pointers to vector data */
138 yy_data = N_VGetArrayPointer(yy);
139 yp_data = N_VGetArrayPointer(yp);
140 rr_data = N_VGetArrayPointer(rr);
141 ewtdata = N_VGetArrayPointer(F2C_IDA_ewtvec);
142 vdata = N_VGetArrayPointer(v);
143 Jvdata = N_VGetArrayPointer(Jv);
144 v1data = N_VGetArrayPointer(vtemp1);
145 v2data = N_VGetArrayPointer(vtemp2);
146
147 IDA_userdata = (FIDAUserData) user_data;
148
149 /* Call user-supplied routine */
150 FIDA_JTIMES(&t, yy_data, yp_data, rr_data, vdata, Jvdata,
151 &c_j, ewtdata, &h,
152 IDA_userdata->ipar, IDA_userdata->rpar,
153 v1data, v2data, &ier);
154
155 return(ier);
156 }
157