/[ascend]/trunk/solvers/ida/idalinear.c
ViewVC logotype

Contents of /trunk/solvers/ida/idalinear.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1507 - (show annotations) (download) (as text)
Wed Jun 27 11:25:37 2007 UTC (15 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 7126 byte(s)
Moving integrators to own directory, about to make them self-contained shared libraries.
1 #include "idalinear.h"
2 #include <ida/ida_impl.h>
3 #include <utilities/error.h>
4 #include <utilities/ascMalloc.h>
5
6 #include <sundials/sundials_math.h>
7 #define ZERO RCONST(0.0)
8 #define ONE RCONST(1.0)
9 #define TWO RCONST(2.0)
10
11 typedef struct IntegratorIdaAscendMemStruct{
12 long integ_neq; /* problem size */
13 IntegratorSparseJacFn *integ_jacfn; /* sparse mtx jacobian evaluation function */
14 void * integ_jac_data; /* data for use by jacobvian evaluation function */
15 int integ_lastflag;
16 unsigned long integ_nje;
17 unsigned long integ_nre;
18 mtx_matrix_t integ_sparse_jac_matrix;
19 } IntegratorIdaAscendMem;
20
21 /* readability replacements (see also ida_dense.c from SUNDIALS distro */
22
23 #define linit (IDA_mem->ida_linit)
24 #define lsetup (IDA_mem->ida_lsetup)
25 #define lsolve (IDA_mem->ida_lsolve)
26 #define lperf (IDA_mem->ida_lperf)
27 #define lfree (IDA_mem->ida_lfree)
28 #define lmem (IDA_mem->ida_lmem)
29 #define tn (IDA_mem->ida_tn)
30 #define cjratio (IDA_mem->ida_cjratio)
31 #define cj (IDA_mem->ida_cj)
32
33 #define nje (iamem->integ_nje)
34 #define nre (iamem->integ_nre)
35 #define lastflag (iamem->integ_lastflag)
36 #define neq (iamem->integ_neq)
37 #define jacfn (iamem->integ_jacfn)
38 #define jacdata (iamem->integ_jac_data)
39 #define JJ (iamem->integ_sparse_jac_matrix)
40
41 #define MSGD_IDAMEM_NULL "Integrator memory is NULL."
42 #define MSGD_MEM_FAIL "A memory request failed."
43 #define MSGD_LMEM_NULL "IDAASCEND memory is NULL."
44 #define MSGD_JACFN_UNDEF "The sparse jacobian evaluation routine has not been provided."
45
46 /*------------------------------------
47 Internal setup/evaluation routines required by IDA. As documented in IDA Manual Ch. 8
48 */
49 int integrator_ida_linit(IDAMem ida_mem);
50
51 int integrator_ida_lsetup(IDAMem ida_mem, N_Vector yyp, N_Vector ypp,
52 N_Vector resp, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3
53 );
54
55 int integrator_ida_lsolve(IDAMem ida_mem, N_Vector b, N_Vector weight,
56 N_Vector ycur, N_Vector ypcur, N_Vector rescur
57 );
58
59 int integrator_ida_lfree(IDAMem ida_mem);
60
61 /*------------------------------------
62 User functions (called from ida.c in this directory)
63 */
64
65 int IDAASCEND(void *ida_mem, long _neq){
66 IDAMem IDA_mem;
67 IntegratorIdaAscendMem *iamem;
68
69 if(ida_mem == NULL){
70 IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
71 return(IDAASCEND_MEM_NULL);
72 }
73 IDA_mem = (IDAMem)ida_mem;
74
75 iamem = ASC_NEW(IntegratorIdaAscendMem);
76 if(iamem == NULL){
77 return IDAASCEND_MEM_FAIL;
78 }
79 lmem = (void *)iamem;
80
81 /* free linsolver memory with the previous lfree fn, if allocated */
82 if(lfree != NULL)lfree(ida_mem);
83
84 /* set the internal-use linear solver function pointers for IDA */
85 linit = &integrator_ida_linit;
86 lsetup = &integrator_ida_lsetup;
87 lsolve = &integrator_ida_lsolve;
88 lperf = NULL;
89 lfree = &integrator_ida_lfree;
90
91 /* no jacobian assigned, initially (we will throw an error if the user doesn't assign it though) */
92 jacfn = NULL;
93 lastflag = IDAASCEND_SUCCESS;
94 neq = _neq;
95
96 /* allocate memory required for Jacobian mtx_matrix_t ?? */
97
98 return IDAASCEND_SUCCESS;
99 }
100
101 int IDAASCENDSetJacFn(void *ida_mem, IntegratorSparseJacFn *_jacfn, void *_jac_data){
102 IDAMem IDA_mem;
103 IntegratorIdaAscendMem *iamem;
104
105 if(ida_mem == NULL){
106 IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
107 return(IDAASCEND_MEM_NULL);
108 }
109 IDA_mem = (IDAMem)ida_mem;
110
111 if(lmem == NULL){
112 IDAProcessError(ida_mem, IDAASCEND_LMEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_LMEM_NULL);
113 return(IDAASCEND_LMEM_NULL);
114 }
115 iamem = (IntegratorIdaAscendMem *)lmem;
116
117 jacfn = _jacfn;
118
119 return IDAASCEND_SUCCESS;
120 }
121
122 int IDAASCENDGetLastFlag(void *ida_mem, int *flag){
123 IDAMem IDA_mem;
124 IntegratorIdaAscendMem *iamem;
125
126 if(ida_mem == NULL){
127 IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
128 return(IDAASCEND_MEM_NULL);
129 }
130 IDA_mem = (IDAMem)ida_mem;
131
132 if(lmem == NULL){
133 IDAProcessError(ida_mem, IDAASCEND_LMEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_LMEM_NULL);
134 return(IDAASCEND_LMEM_NULL);
135 }
136 iamem = (IntegratorIdaAscendMem *)lmem;
137
138 *flag = lastflag;
139
140 return IDAASCEND_SUCCESS;
141 }
142
143 char *IDAASCENDGetReturnFlagName(int flag){
144 char *name;
145
146 name = ASC_NEW_ARRAY(char,30);
147
148 #define FLAG(N) case N:sprintf(name,#N);break
149
150 switch(flag) {
151 FLAG(IDAASCEND_SUCCESS);
152 FLAG(IDAASCEND_MEM_NULL);
153 FLAG(IDAASCEND_LMEM_NULL);
154 FLAG(IDAASCEND_MEM_FAIL);
155 FLAG(IDAASCEND_JACFN_UNDEF);
156 FLAG(IDAASCEND_JACFN_UNRECVR);
157 FLAG(IDAASCEND_JACFN_RECVR);
158 default:
159 sprintf(name,"Unknown flag value '%d'",flag);
160 }
161
162 return name;
163 }
164
165
166 /*------------------------------------
167 Internal setup/evaluation routines required by IDA
168 */
169
170 int integrator_ida_linit(IDAMem IDA_mem){
171 IntegratorIdaAscendMem *iamem;
172 iamem = (IntegratorIdaAscendMem *)lmem;
173
174 CONSOLE_DEBUG("Initialising IDA linear solver");
175 nje = 0;
176 nre = 0;
177 jacfn = NULL;
178
179 /* initialise anything else in the IntegratorIdaAscendMem that needs it */
180
181 return 0;
182 }
183
184 int integrator_ida_lsetup(IDAMem IDA_mem
185 , N_Vector yyp, N_Vector ypp, N_Vector rrp
186 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
187 ){
188 int retval;
189
190 IntegratorIdaAscendMem *iamem;
191 iamem = (IntegratorIdaAscendMem *)lmem;
192
193 CONSOLE_DEBUG("Setting up IDA linear problem");
194
195 if(jacfn==NULL){
196 lastflag = IDAASCEND_JACFN_UNDEF;
197 return -1; /* unrecoverable */
198 }
199
200 /* Increment nje counter. */
201 nje++;
202
203 /* clear the jacobian matrix */
204 /* ... */
205
206 /* evaluate the jacobian */
207 retval = jacfn(neq, tn, yyp, ypp, rrp, cj, jacdata, JJ,
208 tmp1, tmp2, tmp3
209 );
210
211 if(retval < 0){
212 lastflag = IDAASCEND_JACFN_UNRECVR;
213 return -1;
214 }
215 if (retval > 0) {
216 lastflag = IDAASCEND_JACFN_RECVR;
217 return +1;
218 }
219
220 /* do block decomposition, LU factorisation or whatever, return success or fail flag */
221 /* ... */
222
223 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
224 return(-1);
225 }
226
227 /**
228 This routines handles the linear solve operation for the IDAASCEND linear
229 solver. It interfaces to the appropriate mtx_matrix routines for this,
230 but also scales solution vector according to cjratio.
231
232 @return IDAASCEND_SUCESS on success
233 */
234 int integrator_ida_lsolve(IDAMem IDA_mem
235 , N_Vector b, N_Vector weight
236 , N_Vector ycur, N_Vector ypcur, N_Vector rrcur
237 ){
238 realtype *bd;
239
240 IntegratorIdaAscendMem *iamem;
241 iamem = (IntegratorIdaAscendMem *)lmem;
242
243 /* retrieve the data array for the RHS vector, 'b' */
244 bd = N_VGetArrayPointer(b);
245
246 /* call the necessary linsolqr routine here */
247 /* ... */
248 /* For IDADENSE, it's: DenseGETRS(JJ, pivots, bd); */
249
250 /* not sure what this is doing yet */
251 /* For IDADENSE, it's: Scale the correction to account for change in cj. */
252 if(cjratio != ONE){
253 N_VScale(TWO/(ONE + cjratio), b, b);
254 }
255
256 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
257 return -1;
258 }
259
260 int integrator_ida_lfree(IDAMem IDA_mem){
261 CONSOLE_DEBUG("Freeing IDA linear solver data");
262
263 /* free jacobian mtx_matrix_t data ?? */
264
265 if(lmem!=NULL){
266 ASC_FREE(lmem);
267 lmem=NULL;
268 }
269 return 0;
270 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22