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

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