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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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