/[ascend]/trunk/base/generic/solver/lsode.c
ViewVC logotype

Annotation of /trunk/base/generic/solver/lsode.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1133 - (hide annotations) (download) (as text)
Sun Jan 14 11:51:48 2007 UTC (15 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 38489 byte(s)
Removed the 'REX' and 'IEX' array aliases (added the solver parameter comments directly in the slv_param_* calls).
Renamed Simulation::[gs]etSolverParameters to Simulation::[gs]etParameters.
Added [gs]etParameter (singular) methods to SWIG wrapper.
Fixed missing NULL in IDA.
Fixed bug with BDF/AM wrong way round in LSODE.
Removed some debug output from slv.c
1 johnpye 1132 /* :ex: set ts=2 */
2 johnpye 669 /* ASCEND modelling environment
3     Copyright 1997, Carnegie Mellon University
4 johnpye 1129 Copyright (C) 2006-2007 Carnegie Mellon University
5 johnpye 669
6     This program is free software; you can redistribute it and/or modify
7     it under the terms of the GNU General Public License as published by
8     the Free Software Foundation; either version 2, or (at your option)
9     any later version.
10    
11     This program is distributed in the hope that it will be useful,
12     but WITHOUT ANY WARRANTY; without even the implied warranty of
13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14     GNU General Public License for more details.
15    
16     You should have received a copy of the GNU General Public License
17     along with this program; if not, write to the Free Software
18     Foundation, Inc., 59 Temple Place - Suite 330,
19     Boston, MA 02111-1307, USA.
20 johnpye 1129 *//** @file
21     LSODE Integrator
22 johnpye 669
23 johnpye 1129 Here we are solving the system of first order ordinary differential
24     equations, ydot = F(y,t), with y(t_0)=y_0 given.
25    
26     In some places 'x' is named instead of 't'. We're trying to migrate to using
27     't' to label the independent variable.
28    
29     There is some excellent documentation about the LSODE integrator available
30     in the document "Description and Use of LSODE, the Livermore Solver for
31     Ordinary Differential Equations, by K Radhakrishnan and A C Hindmarsh.
32     http://www.llnl.gov/CASC/nsde/pubs/u113855.pdf
33    
34 johnpye 669 (old implementation notes:)
35    
36     As fortran io is unreliably portable (vc5+digital fortran)
37     we have converted xerrwv to xascwv provided here.
38 johnpye 741
39 johnpye 669 The lsode interface variable t is actually an array of
40     2 doubles rather than just 1. The first is the the one
41     used by lsode. The second is used by LSODE_FEX to tell
42     what the last time it was called was. This is so the
43     C driver can tell if it needs to resolve d to compute
44     observation variables. If x[0]==x[1] we save ourselves
45     a solve.
46    
47 johnpye 741 @NOTE The above doesn't work since lsode doesn't use the same t internally
48 johnpye 1129 that we hand it. @ENDNOTE
49 johnpye 669 *//*
50 johnpye 1129 original version by Kirk Abbott and Ben Allan, Jan 1994.
51     Last in CVS $Revision: 1.29 $ $Date: 2000/01/25 02:26:31 $ $Author: ballan $
52     */
53 johnpye 669
54     #ifndef NO_SIGNAL_TRAPS
55     #include <signal.h>
56     #include <setjmp.h>
57     #endif /* NO_SIGNAL_TRAPS */
58    
59     #include <utilities/ascConfig.h>
60     #include <utilities/error.h>
61     #include <compiler/instance_enum.h>
62     #include <utilities/ascSignal.h>
63     #include <utilities/ascMalloc.h>
64 johnpye 888 #include <utilities/ascPanic.h>
65 johnpye 669
66     #include "slv_types.h"
67     #include "mtx.h"
68     #include "rel.h"
69     #include "var.h"
70     #include "discrete.h"
71     #include "conditional.h"
72     #include "bnd.h"
73     #include "logrel.h"
74     #include "slv_common.h"
75     #include "linsol.h"
76     #include "linsolqr.h"
77     #include "slv_client.h"
78 johnpye 1055 #include "slv_stdcalls.h"
79 johnpye 1069 #include <packages/sensitivity.h>
80 johnpye 1128 #include "densemtx.h"
81 johnpye 669
82     #include "integrator.h"
83     #include "lsode.h"
84    
85 johnpye 977 const IntegratorInternals integrator_lsode_internals = {
86     integrator_lsode_create
87     ,integrator_lsode_params_default
88     ,integrator_analyse_ode /* note, this routine is back in integrator.c */
89     ,integrator_lsode_solve
90 johnpye 1129 ,integrator_lsode_write_matrix
91 johnpye 977 ,integrator_lsode_free
92     ,INTEG_LSODE
93     ,"LSODE"
94     };
95    
96 johnpye 669 /*
97     #include "Sensitivity.h"
98     *//* see the packages dir */
99    
100     /*
101     * NOUNDERBARS --> FORTRAN compiler naming convention for subroutine
102     * is wierd. WIN32/CRAY is treated as special case
103     */
104     #ifdef APOLLO
105     #define NOUNDERBARS TRUE
106     #endif
107     #ifdef _HPUX_SOURCE
108     #define NOUNDERBARS TRUE
109     #endif
110     /* AIX xlf will not suffix an underbar on a symbol
111     * unless xlf is given the ``-qextname'' option
112     */
113     #ifdef _AIX
114     #define NOUNDERBARS TRUE
115     #endif
116    
117     #ifdef NOUNDERBARS
118     #define LSODE lsode
119     #define LSODE_JEX jex
120     #define LSODE_FEX fex
121     #define GETCOMMON get_lsode_common
122     #define XASCWV xascwv
123     #else
124     /* sun, __alpha, __sgi, ... */
125     #define LSODE lsode_
126     #define LSODE_JEX jex_
127     #define LSODE_FEX fex_
128     #define GETCOMMON get_lsode_common_
129     #define XASCWV xascwv_
130     #endif
131    
132     #if defined(CRAY) || (defined(__WIN32__) && !defined(__MINGW32_VERSION))
133     #undef LSODE
134     #undef LSODE_JEX
135     #undef LSODE_FEX
136     #undef GETCOMMON
137     #undef XASCWV
138     #define XASCWV XASCWV
139     #define LSODE LSODE
140     #define LSODE_JEX JEX
141     #define LSODE_FEX FEX
142     #define GETCOMMON GET_LSODE_COMMON
143     #endif
144    
145     #define DOTIME FALSE
146    
147     /* definitions of lsode supported children of atoms, etc */
148     /********************************************************************/
149     /* solver_var children expected for state variables */
150     static symchar *g_symbols[2];
151     #define STATERTOL g_symbols[0]
152     #define STATEATOL g_symbols[1]
153     static
154     void InitTolNames(void)
155     {
156     STATERTOL = AddSymbol("ode_rtol");
157     STATEATOL = AddSymbol("ode_atol");
158     }
159    
160 johnpye 1097 /*--------------------------
161     Data space for use by LSODE
162     */
163    
164 johnpye 669 /**
165 johnpye 741 Because LSODE doesn't seem to make an allowance for 'client data' we
166 johnpye 669 have to store this as a 'local global' and fish it out when we're in the
167     callbacks.
168     */
169 johnpye 1097 static IntegratorSystem *l_lsode_blsys = 0;
170 johnpye 669
171 johnpye 1097 typedef enum{
172     lsode_none=0 /* true on first call */
173     , lsode_function /**< a function evaluation */
174     , lsode_derivative /**< a gradient evaluation */
175     } IntegratorLsodeLastCallType;
176     /**
177     Enumeration that tells ASCEND what was the most recent thing LSODE did.
178     */
179 johnpye 669
180 johnpye 1097 typedef enum{
181     lsode_ok=0 /**< success */
182     , lsode_nok /**< bad return from func or grad */
183     } IntegratorLsodeStatusCode;
184     /**
185     Enumeration to tell ASCEND if anything failed in a FEX or JEX call.
186     */
187 johnpye 669
188 johnpye 1128 typedef struct IntegratorLsodeDataStruct{
189 johnpye 1097 long n_eqns; /**< dimension of state vector */
190     int *input_indices; /**< vector of state vars indexes */
191     int *output_indices; /**< vector of derivative var indexes */
192 johnpye 1129 struct var_variable **y_vars; /**< NULL-terminated list of states vars */
193     struct var_variable **ydot_vars; /**< NULL-terminated list of derivative vars*/
194     struct rel_relation **rlist; /**< NULL-terminated list of relevant rels
195 johnpye 1097 to be differentiated */
196 johnpye 1129 DenseMatrix dydot_dy; /**< change in derivatives wrt states */
197 johnpye 1097
198     IntegratorLsodeLastCallType lastcall; /* type of last call; func or grad */
199     IntegratorLsodeStatusCode status; /* solve status */
200     char stop; /* stop requested? */
201     int partitioned; /* partioned func evals or not */
202    
203 johnpye 669 } IntegratorLsodeData;
204    
205 johnpye 1097 /**<
206     Global data structure for LSODE.
207 johnpye 669
208 johnpye 1097 @NOTE LSODE is not reentrant! @ENDNOTE
209     */
210    
211     /** Macro to declare a local var and fetch the 'enginedata' stuff into it from l_lsode_blsys. */
212     #define LSODEDATA_GET(N) \
213     IntegratorLsodeData *N; \
214     asc_assert(l_lsode_blsys!=NULL); \
215     N = (IntegratorLsodeData *)l_lsode_blsys->enginedata; \
216     asc_assert(N!=NULL)
217    
218     /** Macro to set the globa l_lsode_blsys to the currently blsys ptr. */
219     #define LSODEDATA_SET(N) \
220     asc_assert(l_lsode_blsys==NULL); \
221     asc_assert(N!=NULL); \
222     l_lsode_blsys = N
223    
224     #define LSODEDATA_RELEASE() \
225     asc_assert(l_lsode_blsys!=NULL); \
226     l_lsode_blsys = NULL;
227    
228 johnpye 669 /*----------------------------
229     Function types that LSODE wants to use
230     */
231    
232     /**
233     Type of function used to evaluate derivative system.
234     */
235 johnpye 888 typedef void LsodeEvalFn(int *, double *, double *, double *);
236 johnpye 669
237     /**
238     Type of function used to evaluate jacobian system.
239     */
240 johnpye 888 typedef void LsodeJacobianFn(int *, double *, double *, int *, int *, double *, int *);
241 johnpye 669
242     /*----------------------------
243     forward declarations
244     */
245    
246     int integrator_lsode_setup_diffs(IntegratorSystem *blsys);
247    
248     /**
249     void LSODE(&fex, &neq, y, &x, &xend, &itol, reltol, abtol, &itask,
250     &istate, &iopt ,rwork, &lrw, iwork, &liw, &jex, &mf);
251    
252     This is a prototype for the *fortran* LSODE function.
253    
254     No 'extern' here, so we want linker to complain if no static linkage.
255     */
256     void LSODE(LsodeEvalFn*,int *neq ,double *y ,double *x
257     ,double *xend
258     ,int *itol ,double *reltol ,double *abtol
259     ,int *itask ,int *istate ,int *iopt
260     ,double *rwork ,int *lrw
261     ,int *iwork ,int *liw
262     ,LsodeJacobianFn *jex ,int *mf
263     );
264    
265     /*------------------------------------------------------
266     Memory allocation/free
267     */
268    
269     void integrator_lsode_create(IntegratorSystem *blsys){
270     IntegratorLsodeData *d;
271     d = ASC_NEW_CLEAR(IntegratorLsodeData);
272     d->n_eqns=0;
273     d->input_indices=NULL;
274     d->output_indices=NULL;
275     d->y_vars=NULL;
276     d->ydot_vars=NULL;
277     d->rlist=NULL;
278 johnpye 1129 d->dydot_dy=DENSEMATRIX_EMPTY;
279 johnpye 669 blsys->enginedata=(void*)d;
280 johnpye 962 integrator_lsode_params_default(blsys);
281    
282 johnpye 669 }
283    
284     /**
285 johnpye 888 Cleanup the data struct that belongs to LSODE
286 johnpye 669 */
287     void integrator_lsode_free(void *enginedata){
288     IntegratorLsodeData d;
289     d = *((IntegratorLsodeData *)enginedata);
290    
291     if(d.input_indices)ASC_FREE(d.input_indices);
292     d.input_indices = NULL;
293    
294     if(d.output_indices)ASC_FREE(d.output_indices);
295     d.output_indices = NULL;
296    
297     if(d.y_vars)ASC_FREE(d.y_vars);
298     d.y_vars = NULL;
299    
300     if(d.ydot_vars)ASC_FREE(d.ydot_vars);
301     d.ydot_vars = NULL;
302    
303     if(d.rlist)ASC_FREE(d.rlist);
304     d.rlist = NULL;
305    
306 johnpye 1129 densematrix_destroy(d.dydot_dy);
307 johnpye 669
308     d.n_eqns = 0L;
309     }
310    
311 johnpye 942 /*------------------------------------------------------------------------------
312     PARAMETERS
313     */
314    
315     enum ida_parameters{
316 johnpye 1132 LSODE_PARAM_METH
317     ,LSODE_PARAM_MITER
318     ,LSODE_PARAM_MAXORD
319     ,LSODE_PARAM_TIMING
320 johnpye 963 ,LSODE_PARAM_RTOLVECT
321     ,LSODE_PARAM_RTOL
322     ,LSODE_PARAM_ATOLVECT
323     ,LSODE_PARAM_ATOL
324 johnpye 942 ,LSODE_PARAMS_SIZE
325     };
326    
327     /**
328     Here the full set of parameters is defined, along with upper/lower bounds,
329     etc. The values are stuck into the blsys->params structure.
330    
331     @return 0 on success
332     */
333     int integrator_lsode_params_default(IntegratorSystem *blsys){
334    
335     asc_assert(blsys!=NULL);
336 johnpye 962 asc_assert(blsys->engine==INTEG_LSODE);
337 johnpye 942 slv_parameters_t *p;
338     p = &(blsys->params);
339    
340     slv_destroy_parms(p);
341    
342     if(p->parms==NULL){
343     p->parms = ASC_NEW_ARRAY(struct slv_parameter, LSODE_PARAMS_SIZE);
344     if(p->parms==NULL)return -1;
345     p->dynamic_parms = 1;
346     }else{
347     asc_assert(p->num_parms == LSODE_PARAMS_SIZE);
348     }
349    
350     /* reset the number of parameters to zero so that we can check it at the end */
351     p->num_parms = 0;
352    
353 johnpye 1132 slv_param_char(p,LSODE_PARAM_METH
354     ,(SlvParameterInitChar){{"meth"
355     ,"Integration method",1
356     ,"AM=Adams-Moulton method (for non-stiff problems), BDF=Backwards"
357     " Difference Formular (for stiff problems). See 'Description and"
358     " Use of LSODE', section 3.1."
359 johnpye 1133 }, "BDF"}, (char *[]){"AM","BDF",NULL}
360 johnpye 1132 );
361    
362     slv_param_int(p,LSODE_PARAM_MITER
363     ,(SlvParameterInitInt){{"miter"
364     ,"Corrector iteration technique", 1
365     ,"0=Functional iteration, 1=Modified Newton iteration with user-"
366     "supplied analytical Jacobian, 2=Modified Newton iteration with"
367     " internally-generated numerical Jacobian, 3=Modified Jacobi-Newton"
368     " iteration with internally generated numerical Jacobian. See "
369     " 'Description and Use of LSODE', section 3.1. Note that not all"
370     " methods described there are available via ASCEND."
371     }, 1, 0, 3}
372     );
373    
374     slv_param_int(p,LSODE_PARAM_MAXORD
375     ,(SlvParameterInitInt){{"maxord"
376     ,"Maximum method order", 1
377     ,"See 'Description and Use of LSODE', p 92 and p 8. Limits <=12 for BDF"
378     " and <=5 for AM. Higher values are reduced automatically. See notes on"
379     " page 92 regarding oscillatory systems."
380     }, 12, 1, 12}
381     );
382    
383 johnpye 942 slv_param_bool(p,LSODE_PARAM_TIMING
384     ,(SlvParameterInitBool){{"timing"
385 johnpye 1132 ,"Output timing statistics?",1
386     ,"If TRUE, additional timing statistics will be output to the"
387     " console during integration."
388 johnpye 942 }, TRUE}
389     );
390    
391 johnpye 963 slv_param_bool(p,LSODE_PARAM_ATOLVECT
392     ,(SlvParameterInitBool){{"atolvect"
393     ,"Use 'ode_atol' values as specified for each var?",1
394 johnpye 1132 ,"If TRUE, values of 'ode_atol' are taken from your model and used"
395 johnpye 963 " in the integration. If FALSE, a scalar absolute tolerance (atol)"
396     " is shared by all variables."
397     }, TRUE }
398     );
399    
400     slv_param_real(p,LSODE_PARAM_ATOL
401     ,(SlvParameterInitReal){{"atol"
402 johnpye 962 ,"Scalar absolute error tolerance",1
403     ,"Default value of the scalar absolute error tolerance (for cases"
404     " where not specified in oda_atol var property. See 'lsode.f' for"
405     " details"
406     }, 1e-6, DBL_MIN, DBL_MAX }
407     );
408    
409 johnpye 963 slv_param_bool(p,LSODE_PARAM_RTOLVECT
410     ,(SlvParameterInitBool){{"rtolvect"
411     ,"Use 'ode_rtol' values as specified for each var?",1
412     ,"If TRUE, values of 'ode_atol' are taken from your model and used "
413     " in the integration. If FALSE, a scalar absolute tolerance (rtol)"
414     " is shared by all variables."
415     }, TRUE }
416     );
417    
418     slv_param_real(p,LSODE_PARAM_RTOL
419     ,(SlvParameterInitReal){{"rtol"
420 johnpye 962 ,"Scalar relative error tolerance",1
421     ,"Default value of the scalar relative error tolerance (for cases"
422     " where not specified in oda_rtol var property. See 'lsode.f' for"
423     " details"
424     }, 1e-6, DBL_MIN, DBL_MAX }
425     );
426    
427 johnpye 942 asc_assert(p->num_parms == LSODE_PARAMS_SIZE);
428     return 0;
429     }
430    
431 johnpye 669 /*------------------------------------------------------------------------------
432     PROBLEM ANALYSIS
433     */
434    
435 ben.allan 704 /**
436     @TODO needs work. Assumes struct Instance* and struct var_variable*
437     are synonymous, which demonstrates the need for a method to take
438     an instance and ask the solvers for its global or local index
439     if var and inst are decoupled.
440     */
441     int integrator_lsode_setup_diffs(IntegratorSystem *blsys) {
442 johnpye 908 /* long n_eqns; */
443 ben.allan 704 unsigned long nch,i;
444    
445     struct var_variable **vp;
446 johnpye 669 int *ip;
447    
448     IntegratorLsodeData *enginedata;
449 johnpye 1097 asc_assert(blsys!=NULL);
450 ben.allan 704 enginedata = (IntegratorLsodeData *)blsys->enginedata;
451 johnpye 669
452     assert(enginedata->n_eqns==blsys->n_y);
453 ben.allan 704
454 johnpye 1097 /*
455     Put the
456     Let us now process what we consider *inputs* to the problem as
457     far as ASCEND is concerned; i.e. the state vars or the y_vars's
458     if you prefer.
459     */
460     nch = enginedata->n_eqns;
461 johnpye 669
462 johnpye 1097 vp = enginedata->y_vars;
463     ip = enginedata->input_indices;
464     for (i=0;i<nch;i++) {
465     *vp = (struct var_variable *)blsys->y[i];
466     *ip = var_sindex(*vp);
467     vp++;
468     ip++;
469     }
470     *vp = NULL; /* terminate */
471 ben.allan 704
472 johnpye 1097 /*
473     Let us now go for the outputs, ie the derivative terms.
474     */
475     vp = enginedata->ydot_vars;
476     ip = enginedata->output_indices;
477     for (i=0;i<nch;i++) {
478     *vp = (struct var_variable *)blsys->ydot[i];
479     *ip = var_sindex(*vp);
480     vp++; /* dont assume that a var is synonymous with */
481     ip++; /* an Instance; that might/will change soon */
482     }
483     *vp = NULL; /* terminate */
484 ben.allan 704
485 johnpye 1097 return 0;
486 johnpye 669 }
487    
488     /**
489 johnpye 888 allocates, fills, and returns the atol vector based on LSODE
490 johnpye 669
491 johnpye 963 State variables missing child ode_rtol will be defaulted to ATOL
492 johnpye 669 */
493 johnpye 888 static double *lsode_get_atol( IntegratorSystem *blsys) {
494 johnpye 669
495     struct Instance *tol;
496     double *atoli;
497     int i,len;
498 johnpye 963 double atol;
499 johnpye 669
500     len = blsys->n_y;
501 johnpye 963 atoli = ASC_NEW_ARRAY(double, blsys->n_y); /* changed, this was n_y+1 before, dunnowi -- JP */
502 johnpye 669 if (atoli == NULL) {
503     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
504     return atoli;
505     }
506 johnpye 963
507     if(!SLV_PARAM_BOOL(&(blsys->params),LSODE_PARAM_ATOLVECT)){
508     atol = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_ATOL);
509     CONSOLE_DEBUG("Using ATOL = %f for all vars", atol);
510     for(i=0; i<len; ++i){
511     atoli[i] = atol;
512     }
513     }else{
514     InitTolNames();
515     for (i=0; i<len; i++) {
516    
517     tol = ChildByChar(var_instance(blsys->y[i]),STATEATOL);
518     if (tol == NULL || !AtomAssigned(tol) ) {
519     atoli[i] = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_ATOL);
520     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming atol = %3g"
521     "for ode_atol child undefined for state variable %ld."
522     ,atoli[i], blsys->y_id[i]
523     );
524     } else {
525     atoli[i] = RealAtomValue(tol);
526     CONSOLE_DEBUG("Using atol %3g for state variable %ld.",atoli[i], blsys->y_id[i]);
527     }
528 johnpye 669 }
529     }
530     return atoli;
531     }
532    
533     /**
534 johnpye 888 Allocates, fills, and returns the rtol vector based on LSODE
535 johnpye 669
536 johnpye 963 State variables missing child ode_rtol will be defaulted to RTOL
537 johnpye 669 */
538 johnpye 888 static double *lsode_get_rtol( IntegratorSystem *blsys) {
539 johnpye 669
540     struct Instance *tol;
541 johnpye 963 double rtol, *rtoli;
542 johnpye 669 int i,len;
543    
544     len = blsys->n_y;
545     rtoli = ASC_NEW_ARRAY(double, blsys->n_y+1);
546     if (rtoli == NULL) {
547     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
548     return rtoli;
549     }
550 johnpye 963 if(!SLV_PARAM_BOOL(&(blsys->params),LSODE_PARAM_RTOLVECT)){
551     rtol = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_RTOL);
552     CONSOLE_DEBUG("Using RTOL = %f for all vars", rtol);
553     for(i=0; i<len; ++i){
554     rtoli[i] = rtol;
555     }
556     }else{
557     InitTolNames();
558     for (i=0; i<len; i++) {
559     tol = ChildByChar(var_instance(blsys->y[i]),STATERTOL);
560     if (tol == NULL || !AtomAssigned(tol) ) {
561     rtoli[i] = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_RTOL);
562 johnpye 669
563 johnpye 963 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming rtol = %3g"
564     "for ode_rtol child undefined for state variable %ld."
565     ,rtoli[i], blsys->y_id[i]
566     );
567    
568     } else {
569     rtoli[i] = RealAtomValue(tol);
570     }
571 johnpye 669 }
572     }
573 johnpye 963 rtoli[len] = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_RTOL);
574 johnpye 669 return rtoli;
575     }
576    
577     /*
578     Write out a a status message based on the istate parameter.
579     */
580 johnpye 888 static void lsode_write_istate( int istate) {
581 johnpye 1097 switch (istate) {
582     case -1:
583     FPRINTF(ASCERR,"Excess steps taken on this call"
584     " (perhaps wrong MF)."); break;
585     case -2:
586     FPRINTF(ASCERR,"Excess accuracy requested"
587     " (tolerances too small)."); break;
588     case -3:
589     FPRINTF(ASCERR,"Illegal input detected"
590     " (see console)."); break;
591     case -4:
592     FPRINTF(ASCERR,"Repeated error test failures"
593     " (check all inputs)."); break;
594     case -5:
595     FPRINTF(ASCERR,"Repeated convergence failures"
596     " (perhaps bad Jacobian supplied, or wrong choice of MF or tolerances)."); break;
597     case -6:
598     FPRINTF(ASCERR,"Error weight became zero during problem"
599     " (solution component i vanished, and atol or atol(i) = 0)."); break;
600     case -7:
601     FPRINTF(ASCERR,"Interrupted? User cancelled operation?"); break;
602     case -8:
603     FPRINTF(ASCERR,"Error in nonlinear solver"); break;
604     default:
605     FPRINTF(ASCERR,"Unknown 'istate' error code %d from LSODE.",istate);
606     break;
607     }
608 johnpye 669 }
609    
610     /**
611     Free memory allocated for the LSODE, but first check.
612     */
613 johnpye 888 static void lsode_free_mem(double *y, double *reltol, double *abtol, double *rwork,
614 johnpye 669 int *iwork, double *obs, double *dydx)
615     {
616     if (y != NULL) {
617     ascfree((double *)y);
618     }
619     if (reltol != NULL) {
620     ascfree((double *)reltol);
621     }
622     if (abtol != NULL) {
623     ascfree((double *)abtol);
624     }
625     if (rwork != NULL) {
626     ascfree((double *)rwork);
627     }
628     if (iwork != NULL) {
629     ascfree((int *)iwork);
630     }
631     if (obs != NULL) {
632     ascfree((double *)obs);
633     }
634     if (dydx != NULL) {
635     ascfree((double *)dydx);
636     }
637     }
638    
639 johnpye 1097 /**
640     "Temporary" derivative evaluation routine (pre 1995!).
641    
642     Ben says: "The proper permanent fix for lsode is to dump it in favor of
643     cvode or dassl." (so: see ida.c)
644    
645     @return 0 on success
646    
647     @NOTE It is assumed the system has been solved at the current point. @ENDNOTE
648     */
649     int integrator_lsode_derivatives(IntegratorSystem *blsys
650     , int ninputs
651     , int noutputs
652     ){
653 johnpye 669 static int n_calls = 0;
654 johnpye 1097 linsolqr_system_t linsys; /* stuff for the linear system & matrix */
655 johnpye 669 mtx_matrix_t mtx;
656     int32 capacity;
657     real64 *scratch_vector = NULL;
658     int result=0;
659 johnpye 1097 IntegratorLsodeData *enginedata;
660 johnpye 669
661 johnpye 1097 asc_assert(blsys!=NULL);
662     enginedata = (IntegratorLsodeData *)blsys->enginedata;
663     asc_assert(enginedata!=NULL);
664 johnpye 1129 asc_assert(DENSEMATRIX_DATA(enginedata->dydot_dy)!=NULL);
665 johnpye 1097 asc_assert(enginedata->input_indices!=NULL);
666    
667 johnpye 1129 double **dydot_dy = DENSEMATRIX_DATA(enginedata->dydot_dy);
668 johnpye 1097 int *inputs_ndx_list = enginedata->input_indices;
669     int *outputs_ndx_list = enginedata->output_indices;
670     asc_assert(ninputs == blsys->n_y);
671    
672 johnpye 669 (void)NumberFreeVars(NULL); /* used to re-init the system */
673     (void)NumberIncludedRels(NULL); /* used to re-init the system */
674 johnpye 1097 if (!blsys->system) {
675 johnpye 669 FPRINTF(stderr,"The solve system does not exist !\n");
676     return 1;
677     }
678    
679 johnpye 1097 result = Compute_J(blsys->system);
680 johnpye 669 if (result) {
681     FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
682     return 1;
683     }
684    
685 johnpye 1097 linsys = slv_get_linsolqr_sys(blsys->system); /* get the linear system */
686     if (linsys==NULL) {
687 johnpye 669 FPRINTF(stderr,"Early termination due to missing linsolqr system.\n");
688     return 1;
689     }
690 johnpye 1097 mtx = slv_get_sys_mtx(blsys->system); /* get the matrix */
691 johnpye 669 if (mtx==NULL) {
692     FPRINTF(stderr,"Early termination due to missing mtx in linsolqr.\n");
693     return 1;
694     }
695     capacity = mtx_capacity(mtx);
696     scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
697 johnpye 1097 linsolqr_add_rhs(linsys,scratch_vector,FALSE);
698 johnpye 669
699 johnpye 1097 result = LUFactorJacobian(blsys->system);
700 johnpye 669 if (result) {
701     FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
702     goto error;
703     }
704 johnpye 1129 result = Compute_dy_dx_smart(blsys->system, scratch_vector, dydot_dy,
705 johnpye 669 inputs_ndx_list, ninputs,
706     outputs_ndx_list, noutputs);
707    
708 johnpye 1097 linsolqr_remove_rhs(linsys,scratch_vector);
709 johnpye 669 if (result) {
710     FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
711     goto error;
712     }
713    
714 johnpye 1097 error:
715 johnpye 669 n_calls++;
716     if (scratch_vector) {
717     ascfree((char *)scratch_vector);
718     }
719     return result;
720     }
721    
722     /**
723     The current way that we are getting the derivatives (if the problem
724     was solved partitioned) messes up the slv_system so that we *have*
725     to do a *presolve* rather than a simply a *resolve* before doing
726     function calls. This code below attempts to handle these cases.
727     */
728 johnpye 1097 static void LSODE_FEX( int *n_eq ,double *t ,double *y ,double *ydot){
729 johnpye 669 slv_status_t status;
730 johnpye 1097 LSODEDATA_GET(lsodedata);
731 johnpye 669
732     /* slv_parameters_t parameters; pity lsode doesn't allow error returns */
733 johnpye 908 /* int i; */
734 johnpye 1049 unsigned long res;
735 johnpye 669
736     #if DOTIME
737     double time1,time2;
738     #endif
739    
740 johnpye 938 /* CONSOLE_DEBUG("Calling for a function evaluation"); */
741 johnpye 903
742 johnpye 669 #if DOTIME
743     CONSOLE_DEBUG("Calling for a function evaluation");
744     time1 = tm_cpu_time();
745     #endif
746    
747     /*
748 johnpye 888 t[1]=t[0]; can't do this. lsode calls us with a different t than the t we sent in.
749 johnpye 669 */
750     integrator_set_t(l_lsode_blsys, t[0]);
751     integrator_set_y(l_lsode_blsys, y);
752    
753     #if DOTIME
754     time2 = tm_cpu_time();
755     #endif
756    
757 johnpye 1097 switch(lsodedata->lastcall) {
758 johnpye 669 case lsode_none: /* first call */
759 johnpye 903 CONSOLE_DEBUG("FIRST CALL...");
760    
761 johnpye 669 case lsode_derivative:
762 johnpye 1097 if (lsodedata->partitioned) {
763 johnpye 938 /* CONSOLE_DEBUG("PRE-SOLVE"); */
764 johnpye 669 slv_presolve(l_lsode_blsys->system);
765     } else {
766 johnpye 938 /** @TODO this doesn't ever seem to be called */
767 johnpye 903 CONSOLE_DEBUG("RE-SOLVE");
768 johnpye 669 slv_resolve(l_lsode_blsys->system);
769     }
770     break;
771     default:
772     case lsode_function:
773     slv_resolve(l_lsode_blsys->system);
774     break;
775     }
776 johnpye 903
777 johnpye 669 slv_solve(l_lsode_blsys->system);
778     slv_get_status(l_lsode_blsys->system, &status);
779 johnpye 1055 if(slv_check_bounds(l_lsode_blsys->system,0,-1,"")){
780 johnpye 1097 lsodedata->status = lsode_nok;
781 johnpye 1055 }
782    
783 johnpye 888 /* pass the solver status to the integrator */
784 johnpye 1049 res = integrator_checkstatus(status);
785 johnpye 669
786     #if DOTIME
787     time2 = tm_cpu_time() - time2;
788     #endif
789    
790 johnpye 1049 if(res){
791     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for derivatives (%d)",res);
792 johnpye 1055 #if 0
793 johnpye 669 ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
794     FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");
795     for (i = 0; i< *n_eq; i++) {
796     FPRINTF(ASCERR,"y[%4d] = %f\n",i, y[i]);
797     }
798     error_reporter_end_flush();
799 johnpye 1055 #endif
800 johnpye 1097 lsodedata->stop = 1;
801     lsodedata->status = lsode_nok;
802 johnpye 1049 }else{
803 johnpye 1097 lsodedata->status = lsode_ok;
804     /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"lsodedata->status = %d",lsodedata->status); */
805 johnpye 669 }
806     integrator_get_ydot(l_lsode_blsys, ydot);
807    
808 johnpye 1097 lsodedata->lastcall = lsode_function;
809 johnpye 669 #if DOTIME
810     time1 = tm_cpu_time() - time1;
811     CONSOLE_DEBUG("Function evalulation has been completed in time %g. True function call time = %g",time1,time2);
812     #endif
813     }
814    
815     /**
816     Evaluate the jacobian
817     */
818 johnpye 1097 static void LSODE_JEX(int *neq ,double *t, double *y
819     , int *ml ,int *mu ,double *pd, int *nrpd
820     ){
821 johnpye 669 int nok = 0;
822     int i,j;
823    
824 johnpye 1097 LSODEDATA_GET(lsodedata);
825 johnpye 669
826     UNUSED_PARAMETER(t);
827     UNUSED_PARAMETER(y);
828     UNUSED_PARAMETER(ml);
829     UNUSED_PARAMETER(mu);
830    
831 johnpye 938 /* CONSOLE_DEBUG("Calling for a gradient evaluation"); */
832 johnpye 669 #if DOTIME
833     double time1;
834    
835     CONSOLE_DEBUG("Calling for a gradient evaluation");
836     time1 = tm_cpu_time();
837     #endif
838     /*
839     * Make the real call.
840     */
841 johnpye 1097
842     nok = integrator_lsode_derivatives(l_lsode_blsys
843 johnpye 669 , *neq
844     , *nrpd
845     );
846    
847 johnpye 1097 if(nok){
848 johnpye 669 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in computing the derivatives for the system. Failing...");
849 johnpye 1097 lsodedata->status = lsode_nok;
850     lsodedata->lastcall = lsode_derivative;
851     lsodedata->stop = 1;
852 johnpye 669 return;
853 johnpye 1097 }else{
854     lsodedata->status = lsode_ok;
855     lsodedata->lastcall = lsode_derivative;
856 johnpye 669 }
857     /*
858 johnpye 903 Map data from C based matrix to Fortan matrix.
859     We will send in a column major ordering vector for pd.
860     */
861 johnpye 1129 asc_assert(*neq == DENSEMATRIX_NCOLS(lsodedata->dydot_dy));
862     asc_assert(*nrpd == DENSEMATRIX_NROWS(lsodedata->dydot_dy));
863 johnpye 888 for (j=0;j<*neq;j++) { /* loop through columnns */
864     for (i=0;i<*nrpd;i++){ /* loop through rows */
865 johnpye 1129 /* CONSOLE_DEBUG("JAC[r=%d,c=%d]=%f",i,j,lsodedata.dydot_dy[i][j]); */
866     *pd++ = DENSEMATRIX_ELEM(lsodedata->dydot_dy,i,j);
867 johnpye 669 }
868     }
869    
870     #if DOTIME
871     time1 = tm_cpu_time() - time1;
872     CONSOLE_DEBUG("Time to do gradient evaluation %g",time1);
873     #endif
874    
875     return;
876     }
877    
878     /**
879     The public function: here we do the actual integration, I guess.
880 johnpye 966
881 johnpye 1049 Return 0 on success
882 johnpye 669 */
883     int integrator_lsode_solve(IntegratorSystem *blsys
884     , unsigned long start_index, unsigned long finish_index
885     ){
886     slv_status_t status;
887     slv_parameters_t params;
888 johnpye 741 IntegratorLsodeData *d;
889 johnpye 669
890     double x[2];
891     double xend,xprev;
892     unsigned long nsamples, neq;
893     long nobs;
894     int itol, itask, mf, lrw, liw;
895     unsigned long index;
896     int istate, iopt;
897     double * rwork;
898     int * iwork;
899     double *y, *abtol, *reltol, *obs, *dydx;
900     int my_neq;
901 johnpye 902 int reporterstatus;
902 johnpye 1132 const char *method; /* Table 3.1 in D&UoLSODE */
903     int miter; /* Table 3.2 in D&UoLSODE */
904     int maxord; /* page 92 in D&UoLSODE */
905 johnpye 669
906     d = (IntegratorLsodeData *)(blsys->enginedata);
907    
908 ben.allan 704 /* the numer of equations must be equal to blsys->n_y, the number of states */
909     d->n_eqns = blsys->n_y;
910     assert(d->n_eqns>0);
911 johnpye 669
912 ben.allan 704 d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
913     d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
914 johnpye 1129 d->dydot_dy = densematrix_create(d->n_eqns,d->n_eqns);
915 johnpye 669
916     d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);
917     d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);
918    
919     integrator_lsode_setup_diffs(blsys);
920    
921 johnpye 1132 /* LSODE should be OK to deal with any linsol/linsolqr-based solver. But for
922     the moment we restrict to just QRSlv. */
923     if(strcmp(slv_solver_name(slv_get_selected_solver(blsys->system)),"QRSlv") != 0) {
924     ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv must be selected before integration.");
925     return 1;
926     }
927 johnpye 669
928 johnpye 1132 CONSOLE_DEBUG("Solver selected is '%s'",slv_solver_name(slv_get_selected_solver(blsys->system)));
929 johnpye 669
930 johnpye 1132 slv_get_status(blsys->system, &status);
931 johnpye 1097
932 johnpye 669 if (status.struct_singular) {
933     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system is structurally singular.");
934 johnpye 1097 d->status = lsode_nok;
935 johnpye 1049 return 2;
936 johnpye 669 }
937    
938     #if defined(STATIC_LSOD) || defined (DYNAMIC_LSOD)
939    
940     /* here we assume integrators.c is in charge of dynamic loading */
941    
942     slv_get_parameters(blsys->system,&params);
943 johnpye 1097 d->partitioned = 1;
944     d->stop = 0; /* clear 'stop' flag */
945 johnpye 669
946 johnpye 1132 /* read METH and MITER parameters, create MF value */
947     method = SLV_PARAM_CHAR(&(blsys->params),LSODE_PARAM_METH);
948     miter = SLV_PARAM_INT(&(blsys->params),LSODE_PARAM_MITER);
949     maxord = SLV_PARAM_INT(&(blsys->params),LSODE_PARAM_MITER);
950     if(miter < 0 || miter > 3){
951     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unacceptable value '%d' of parameter 'miter'",miter);
952     return 5;
953     }
954     if(strcmp(method,"BDF")==0){
955 johnpye 1133 CONSOLE_DEBUG("method = BDF");
956     mf = 20 + miter;
957 johnpye 1132 if(maxord > 5){
958     maxord = 5;
959     CONSOLE_DEBUG("MAXORD reduced to 5 for BDF");
960     }
961     }else if(strcmp(method,"AM")==0){
962 johnpye 1133 CONSOLE_DEBUG("method = AM");
963     if(maxord > 12){
964     maxord = 12;
965     CONSOLE_DEBUG("MAXORD reduced to 12 for AM");
966     }
967     mf = 10 + miter;
968 johnpye 1132 }else{
969     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unacceptable value '%d' of parameter 'meth'",method);
970     return 5;
971     }
972    
973     CONSOLE_DEBUG("MF = %d",mf);
974    
975 johnpye 669 nsamples = integrator_getnsamples(blsys);
976     if (nsamples <2) {
977     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined.");
978 johnpye 1097 d->status = lsode_nok;
979 johnpye 1049 return 3;
980 johnpye 669 }
981     neq = blsys->n_y;
982     nobs = blsys->n_obs;
983    
984 johnpye 966 /* samplelist_debug(blsys->samples); */
985    
986     /* x[0] = integrator_get_t(blsys); */
987     x[0] = integrator_getsample(blsys, 0);
988 johnpye 669 x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */
989 johnpye 1132
990     switch(mf){
991     case 10: case 20:
992     lrw = 20 + neq * (maxord + 1) + 3 * neq;
993     break;
994     case 11: case 12: case 21: case 22:
995     lrw = 22 + 9*neq + neq*neq;
996     break;
997     case 13: case 23:
998     lrw = 22 + neq * (maxord + 1) + 4 * neq;
999     break;
1000     default:
1001     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unknown size requirements for this value of 'mf'");
1002     return 4;
1003     }
1004    
1005 johnpye 669 rwork = ASC_NEW_ARRAY_CLEAR(double, lrw+1);
1006     liw = 20 + neq;
1007     iwork = ASC_NEW_ARRAY_CLEAR(int, liw+1);
1008     y = integrator_get_y(blsys, NULL);
1009 johnpye 888 reltol = lsode_get_rtol(blsys);
1010     abtol = lsode_get_atol(blsys);
1011 johnpye 669 obs = integrator_get_observations(blsys, NULL);
1012     dydx = ASC_NEW_ARRAY_CLEAR(double, neq+1);
1013     if (!y || !obs || !abtol || !reltol || !rwork || !iwork || !dydx) {
1014 johnpye 888 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1015     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for lsode.");
1016 johnpye 1097 d->status = lsode_nok;
1017 johnpye 1049 return 4;
1018 johnpye 669 }
1019    
1020     /*
1021     Prepare args and call lsode.
1022     */
1023     itol = 4;
1024     itask = 1;
1025     istate = 1;
1026     iopt = 1;
1027     rwork[4] = integrator_get_stepzero(blsys);
1028     rwork[5] = integrator_get_maxstep(blsys);
1029     rwork[6] = integrator_get_minstep(blsys);
1030     iwork[5] = integrator_get_maxsubsteps(blsys);
1031 johnpye 1132 iwork[4] = maxord;
1032 johnpye 669
1033 johnpye 966 if(x[0] > integrator_getsample(blsys, 2)){
1034     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value");
1035 johnpye 1049 return 5;
1036 johnpye 966 }
1037    
1038 johnpye 669 /* put the values from derivative system into the record */
1039     integrator_setsample(blsys, start_index, x[0]);
1040    
1041     integrator_output_init(blsys);
1042    
1043 johnpye 979 /* -- store the initial values of all the stuff */
1044     integrator_output_write(blsys);
1045     integrator_output_write_obs(blsys);
1046    
1047 johnpye 669 my_neq = (int)neq;
1048    
1049     /*
1050     First time entering lsode, x is input. After that,
1051     lsode uses x as output (y output is y(x)). To drive
1052     the loop ahead in time, all we need to do is keep upping
1053     xend.
1054     */
1055 johnpye 966
1056 johnpye 669 blsys->currentstep = 0;
1057     for (index = start_index; index < finish_index; index++, blsys->currentstep++) {
1058     xend = integrator_getsample(blsys, index+1);
1059     xprev = x[0];
1060 johnpye 966 asc_assert(xend > xprev);
1061     /* CONSOLE_DEBUG("LSODE call #%lu: x = [%f,%f]", index,xprev,xend); */
1062 johnpye 669
1063     # ifndef NO_SIGNAL_TRAPS
1064 johnpye 997 if (SETJMP(g_fpe_env)==0) {
1065 johnpye 669 # endif /* NO_SIGNAL_TRAPS */
1066    
1067 johnpye 895 /* CONSOLE_DEBUG("Calling LSODE with end-time = %f",xend); */
1068     /*
1069     switch(mf){
1070 johnpye 888 case 10:
1071     CONSOLE_DEBUG("Non-stiff (Adams) method; no Jacobian will be used"); break;
1072     case 21:
1073     CONSOLE_DEBUG("Stiff (BDF) method, user-supplied full Jacobian"); break;
1074     case 22:
1075     CONSOLE_DEBUG("Stiff (BDF) method, internally generated full Jacobian"); break;
1076     case 24:
1077     CONSOLE_DEBUG("Stiff (BDF) method, user-supplied banded jacobian"); break;
1078     case 25:
1079     CONSOLE_DEBUG("Stiff (BDF) method, internally generated banded jacobian"); break;
1080     default:
1081     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid method id %d for LSODE",mf);
1082 johnpye 895 return 0; * failure *
1083 johnpye 888 }
1084 johnpye 895 */
1085 johnpye 669
1086 johnpye 1097 /* provides some rudimentary locking to prevent reentrance*/
1087     LSODEDATA_SET(blsys);
1088    
1089 johnpye 669 LSODE(&(LSODE_FEX), &my_neq, y, x, &xend,
1090     &itol, reltol, abtol, &itask, &istate,
1091     &iopt ,rwork, &lrw, iwork, &liw, &(LSODE_JEX), &mf);
1092    
1093 johnpye 1097 /* clear the global var */
1094     LSODEDATA_RELEASE();
1095 johnpye 669
1096     # ifndef NO_SIGNAL_TRAPS
1097 johnpye 1097 }else{
1098 johnpye 1099 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in LSODE call.");
1099 johnpye 888 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1100 johnpye 1097 d->status = lsode_ok; /* clean up before we go */
1101     d->lastcall = lsode_none;
1102 johnpye 1049 return 6;
1103 johnpye 669 }
1104     # endif /* NO_SIGNAL_TRAPS */
1105    
1106     /* CONSOLE_DEBUG("AFTER %lu LSODE CALL\n", index); */
1107     /* this check is better done in fex,jex, but lsode takes no status */
1108 johnpye 1097 /* if (Solv_C_CheckHalt()) {
1109 johnpye 669 if (istate >= 0) {
1110     istate=-7;
1111     }
1112     }
1113 johnpye 1097 */
1114     if(d->stop){
1115     istate=-8;
1116     }
1117 johnpye 669
1118     if (istate < 0 ) {
1119     /* some kind of error occurred... */
1120     ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
1121 johnpye 888 lsode_write_istate(istate);
1122     FPRINTF(ASCERR, "\nFurthest point reached was t = %g.\n",x[0]);
1123 johnpye 669 error_reporter_end_flush();
1124    
1125 johnpye 888 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1126 johnpye 669 integrator_output_close(blsys);
1127 johnpye 1049 return 7;
1128 johnpye 669 }
1129    
1130 johnpye 1097 if (d->status==lsode_nok) {
1131 johnpye 669 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to an error in derivative computations.");
1132 johnpye 888 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1133 johnpye 1097 d->status = lsode_ok; /* clean up before we go */
1134     d->lastcall = lsode_none;
1135 johnpye 669 integrator_output_close(blsys);
1136 johnpye 1049 return 8;
1137 johnpye 669 }
1138    
1139     integrator_setsample(blsys, index+1, x[0]);
1140     /* record when lsode actually came back */
1141     integrator_set_t(blsys, x[0]);
1142     integrator_set_y(blsys, y);
1143     /* put x,y in d in case lsode got x,y by interpolation, as it does */
1144    
1145 johnpye 902 reporterstatus = integrator_output_write(blsys);
1146 johnpye 669
1147 johnpye 902 if(reporterstatus==0){
1148     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration cancelled");
1149     lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1150 johnpye 1097 d->status = lsode_ok;
1151     d->lastcall = lsode_none;
1152 johnpye 902 integrator_output_close(blsys);
1153 johnpye 1049 return 9;
1154 johnpye 902 }
1155    
1156 johnpye 669 if (nobs > 0) {
1157     # ifndef NO_SIGNAL_TRAPS
1158 johnpye 997 if (SETJMP(g_fpe_env)==0) {
1159 johnpye 669 # endif /* NO_SIGNAL_TRAPS */
1160    
1161     /* solve for obs since d isn't necessarily already
1162     computed there though lsode's x and y may be.
1163     Note that since lsode usually steps beyond xend
1164     x1 usually wouldn't be x0 precisely if the x1/x0
1165     scheme worked, which it doesn't anyway. */
1166    
1167 johnpye 1097 LSODEDATA_SET(blsys);
1168 johnpye 669 LSODE_FEX(&my_neq, x, y, dydx);
1169 johnpye 1097 LSODEDATA_RELEASE();
1170 johnpye 669
1171     /* calculate observations, if any, at returned x and y. */
1172     obs = integrator_get_observations(blsys, obs);
1173    
1174     integrator_output_write_obs(blsys);
1175    
1176     # ifndef NO_SIGNAL_TRAPS
1177     } else {
1178     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in LSODE FEX call.");
1179 johnpye 888 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1180 johnpye 1097 d->status = lsode_ok; /* clean up before we go */
1181     d->lastcall = lsode_none;
1182 johnpye 669 integrator_output_close(blsys);
1183 johnpye 1049 return 10;
1184 johnpye 669 }
1185     # endif /* NO_SIGNAL_TRAPS */
1186     }
1187 johnpye 855 /* CONSOLE_DEBUG("Integration completed from %3g to %3g.",xprev,x[0]); */
1188 johnpye 669 }
1189    
1190 johnpye 888 CONSOLE_DEBUG("...");
1191 johnpye 669 CONSOLE_DEBUG("Number of steps taken: %1d.", iwork[10]);
1192     CONSOLE_DEBUG("Number of function evaluations: %1d.", iwork[11]);
1193     CONSOLE_DEBUG("Number of Jacobian evaluations: %1d.", iwork[12]);
1194 johnpye 888 CONSOLE_DEBUG("...");
1195 johnpye 669
1196    
1197 johnpye 888 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1198 johnpye 669
1199     /*
1200     * return the system to its original state.
1201     */
1202    
1203 johnpye 1097 d->status = lsode_ok;
1204     d->lastcall = lsode_none;
1205 johnpye 669
1206     integrator_output_close(blsys);
1207    
1208 johnpye 888 CONSOLE_DEBUG("--- LSODE done ---");
1209 johnpye 1049 return 0; /* success */
1210 johnpye 669
1211     #else /* STATIC_LSOD || DYNAMIC_LSOD */
1212    
1213     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration will not be performed. LSODE binary not available.");
1214 johnpye 1097 d->status = lsode_nok;
1215 johnpye 1049 return 11;
1216 johnpye 669
1217     #endif
1218     }
1219    
1220     /**
1221 johnpye 888 Function XASCWV is an error reporting function replacing the XERRWV
1222     routine in lsode.f. The call signature is the same with the original Fortran
1223     function.
1224 johnpye 669
1225 johnpye 888 @see the comments for 'xerrwv' from lsode.f, with which XASCWV is compatible...
1226 johnpye 669
1227 johnpye 888 @param msg = the message (hollerith literal or integer array).
1228     @param nmes = the length of msg (number of characters).
1229     @param nerr = the error number (not used).
1230     @param level = the error level..
1231 johnpye 669 0 or 1 means recoverable (control returns to caller).
1232     2 means fatal (run is aborted--see note below).
1233 johnpye 888 @param ni = number of integers (0, 1, or 2) to be printed with message.
1234     @param i1,i2 = integers to be printed, depending on ni.
1235     @param nr = number of reals (0, 1, or 2) to be printed with message.
1236     @param r1,r2 = reals to be printed, depending on nr.
1237     */
1238     void XASCWV( char *msg, /* pointer to start of message */
1239     int *nmes, /* the length of msg (number of characters) */
1240     int *nerr, /* the error number (not used). */
1241 johnpye 908 int *level,
1242 johnpye 669 int *ni,
1243     int *i1,
1244     int *i2,
1245     int *nr,
1246     double *r1,
1247     double *r2
1248 johnpye 888 ){
1249 johnpye 966 static double r1last;
1250    
1251     asc_assert(*level!=2); // LSODE doesn't give level 2 in our version.
1252    
1253 johnpye 888 switch(*nerr){
1254 johnpye 1132 case 17:
1255     if(*ni==2){
1256     ERROR_REPORTER_HERE(ASC_PROG_ERR,"rwork length needed, lenrw = %d > %d = lrw",*i1, *i2);
1257     return;
1258     } break;
1259 johnpye 966 case 52:
1260     if(*nr==2){
1261     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Illegal t = %f, not in range (t - hu,t) = (%f,%f)", r1last, *r1, *r2);
1262     return;
1263     }else if(*nr==1){
1264     r1last = *r1;
1265     return;
1266     } break;
1267 johnpye 1132 case 201:
1268     if(*nr==0 && *ni==0)return;
1269     if(*nr==1 && *ni==1){
1270     ERROR_REPORTER_HERE(ASC_PROG_ERR,"At current t=%f, mxstep=%f steps taken on this call before reaching tout.",*r1,*i1);
1271     return;
1272     } break;
1273 johnpye 888 case 204:
1274 johnpye 1049 if(*nr==0 && *ni==0)return;
1275 johnpye 966 if(*nr==2){
1276     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);
1277     return;
1278     } break;
1279 johnpye 888 case 205:
1280 johnpye 1132 if(*nr==0 && *ni==0)return;
1281 johnpye 966 if(*nr==2){
1282     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Corrector convergence test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);
1283     return;
1284     } break;
1285     case 27:
1286     if(*nr==1 && *ni==1){
1287     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Trouble with INTDY: itask = %d, tout = %f", *i1, *r1);
1288     return;
1289     } break;
1290     }
1291 johnpye 908
1292 johnpye 966 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1293 johnpye 669
1294 johnpye 966 /* note that %.*s means that a string length (integer) and string pointer are being required */
1295     FPRINTF(stderr,"LSODE error: (%d) %.*s",*nerr,*nmes,msg);
1296 johnpye 1132 if(*ni == 1) {
1297     FPRINTF(stderr,"\nwhere i1 = %d",*i1);
1298     }else if(*ni == 2) {
1299     FPRINTF(stderr,"\nwhere i1 = %d, i2 = %d",*i1,*i2);
1300 johnpye 888 }
1301 johnpye 1132 if(*nr == 1) {
1302     FPRINTF(stderr,"\nwhere r1 = %.13g", *r1);
1303     }else if(*nr == 2) {
1304     FPRINTF(stderr,"\nwhere r1 = %.13g, r2 = %.13g", *r1,*r2);
1305 johnpye 888 }
1306     error_reporter_end_flush();
1307 johnpye 669 }
1308 johnpye 1129
1309     int integrator_lsode_write_matrix(const IntegratorSystem *blsys, FILE *fp){
1310     IntegratorLsodeData *enginedata;
1311     asc_assert(blsys!=NULL);
1312     asc_assert(blsys->engine==INTEG_LSODE);
1313     asc_assert(blsys->enginedata);
1314     enginedata = (IntegratorLsodeData *)blsys->enginedata;
1315    
1316     if(!DENSEMATRIX_DATA(enginedata->dydot_dy)){
1317     ERROR_REPORTER_HERE(ASC_PROG_ERR,"dydot_dy contains no data");
1318     }
1319    
1320 johnpye 1131 #ifdef ASC_WITH_MMIO
1321 johnpye 1129 densematrix_write_mmio(enginedata->dydot_dy,fp);
1322     CONSOLE_DEBUG("Returning after matrix output");
1323     return 0;
1324 johnpye 1131 #else
1325     ERROR_REPORTER_HERE(ASC_PROG_ERR,"MMIO routines not available");
1326     return 1;
1327     #endif
1328 johnpye 1129 }
1329    
1330    
1331    
1332    
1333    

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