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

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