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

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