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

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