/[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 1050 - (show annotations) (download) (as text)
Sat Jan 6 04:04:47 2007 UTC (15 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 33861 byte(s)
Adding integration test in TestSteam.testdsgsat.
Naming couple of anonymous structs.
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 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 /*
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 Because LSODE doesn't seem to make an allowance for 'client data' we
152 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 typedef struct IntegratorLsodeDataStruct{
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 } 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 typedef void LsodeEvalFn(int *, double *, double *, double *);
196
197 /**
198 Type of function used to evaluate jacobian system.
199 */
200 typedef void LsodeJacobianFn(int *, double *, double *, int *, int *, double *, int *);
201
202 /*----------------------------
203 forward declarations
204 */
205
206 int integrator_lsode_setup_diffs(IntegratorSystem *blsys);
207 static double **lsode_densematrix_create(int nrows, int ncols);
208 static void lsode_densematrix_destroy(double **matrix,int nrows);
209
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 integrator_lsode_params_default(blsys);
243
244 }
245
246 /**
247 Cleanup the data struct that belongs to LSODE
248 */
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 if(d.dydx_dx)lsode_densematrix_destroy(d.dydx_dx, d.n_eqns);
269 d.dydx_dx = NULL;
270
271 d.n_eqns = 0L;
272 }
273
274 /*------------------------------------------------------------------------------
275 PARAMETERS
276 */
277
278 enum ida_parameters{
279 LSODE_PARAM_TIMING
280 ,LSODE_PARAM_RTOLVECT
281 ,LSODE_PARAM_RTOL
282 ,LSODE_PARAM_ATOLVECT
283 ,LSODE_PARAM_ATOL
284 ,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 asc_assert(blsys->engine==INTEG_LSODE);
297 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 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 ,"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 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 ,"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 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 /*---------------------------------------------------------
365 Couple of matrix methods...?
366 */
367
368 static double **lsode_densematrix_create(int nrows, int ncols){
369 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 static void lsode_densematrix_destroy(double **matrix,int nrows){
381 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 }
391
392 /*------------------------------------------------------------------------------
393 PROBLEM ANALYSIS
394 */
395
396 /**
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 /* long n_eqns; */
404 unsigned long nch,i;
405
406 struct var_variable **vp;
407 int *ip;
408
409 IntegratorLsodeData *enginedata;
410 enginedata = (IntegratorLsodeData *)blsys->enginedata;
411 assert(enginedata!=NULL);
412
413 assert(enginedata->n_eqns==blsys->n_y);
414
415 /*
416 Put the
417 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 nch = enginedata->n_eqns;
422
423
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 }
449
450 /**
451 allocates, fills, and returns the atol vector based on LSODE
452
453 State variables missing child ode_rtol will be defaulted to ATOL
454 */
455 static double *lsode_get_atol( IntegratorSystem *blsys) {
456
457 struct Instance *tol;
458 double *atoli;
459 int i,len;
460 double atol;
461
462 len = blsys->n_y;
463 atoli = ASC_NEW_ARRAY(double, blsys->n_y); /* changed, this was n_y+1 before, dunnowi -- JP */
464 if (atoli == NULL) {
465 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
466 return atoli;
467 }
468
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 }
491 }
492 return atoli;
493 }
494
495 /**
496 Allocates, fills, and returns the rtol vector based on LSODE
497
498 State variables missing child ode_rtol will be defaulted to RTOL
499 */
500 static double *lsode_get_rtol( IntegratorSystem *blsys) {
501
502 struct Instance *tol;
503 double rtol, *rtoli;
504 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 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
525 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 }
534 }
535 rtoli[len] = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_RTOL);
536 return rtoli;
537 }
538
539 /*
540 Write out a a status message based on the istate parameter.
541 */
542 static void lsode_write_istate( int istate) {
543 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 FPRINTF(ASCERR,"Illegal input detected (see console).");
552 break;
553 case -4:
554 FPRINTF(ASCERR,"Repeated error test failures (check all inputs).");
555 break;
556 case -5:
557 FPRINTF(ASCERR,"Repeated convergence failures (perhaps bad Jacobian supplied, or wrong choice of MF or tolerances).");
558 break;
559 case -6:
560 FPRINTF(ASCERR,"Error weight became zero during problem (solution component i vanished, and atol or atol(i) = 0).");
561 break;
562 case -7:
563 FPRINTF(ASCERR,"Interrupted? User cancelled operation?");
564 break;
565 default:
566 FPRINTF(ASCERR,"Unknown 'istate' error code %d from LSODE.",istate);
567 break;
568 }
569 }
570
571 /**
572 Free memory allocated for the LSODE, but first check.
573 */
574 static void lsode_free_mem(double *y, double *reltol, double *abtol, double *rwork,
575 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 * Extended 7/95 baa to deal with linsolqr and lsode.
607 * It is assumed the system has been solved at the current point.
608 *********************************************************************
609 */
610 int lsode_derivatives(slv_system_t sys, double **dy_dx,
611 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 /* int i; */
683 unsigned long res;
684
685 #if DOTIME
686 double time1,time2;
687 #endif
688
689 /* CONSOLE_DEBUG("Calling for a function evaluation"); */
690
691 #if DOTIME
692 CONSOLE_DEBUG("Calling for a function evaluation");
693 time1 = tm_cpu_time();
694 #endif
695
696 /*
697 t[1]=t[0]; can't do this. lsode calls us with a different t than the t we sent in.
698 */
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 CONSOLE_DEBUG("FIRST CALL...");
709
710 case lsode_derivative:
711 if (lsodesys.partitioned) {
712 /* CONSOLE_DEBUG("PRE-SOLVE"); */
713 slv_presolve(l_lsode_blsys->system);
714 } else {
715 /** @TODO this doesn't ever seem to be called */
716 CONSOLE_DEBUG("RE-SOLVE");
717 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
726 slv_solve(l_lsode_blsys->system);
727 slv_get_status(l_lsode_blsys->system, &status);
728 /* pass the solver status to the integrator */
729 res = integrator_checkstatus(status);
730
731 #if DOTIME
732 time2 = tm_cpu_time() - time2;
733 #endif
734
735 if(res){
736 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for derivatives (%d)",res);
737 /*
738 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 */
745 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 /* CONSOLE_DEBUG("Calling for a gradient evaluation"); */
775 #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 nok = lsode_derivatives(l_lsode_blsys->system
785 , 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 Map data from C based matrix to Fortan matrix.
803 We will send in a column major ordering vector for pd.
804 */
805 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 *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
823 Return 0 on success
824 */
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 IntegratorLsodeData *d;
831
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 int reporterstatus;
846
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 /* 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
856 d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
857 d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
858 d->dydx_dx = lsode_densematrix_create(d->n_eqns,d->n_eqns);
859
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 /* 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 1;
870 }
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 2;
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 3;
892 }
893 neq = blsys->n_y;
894 nobs = blsys->n_obs;
895
896 /* samplelist_debug(blsys->samples); */
897
898 /* x[0] = integrator_get_t(blsys); */
899 x[0] = integrator_getsample(blsys, 0);
900 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 reltol = lsode_get_rtol(blsys);
907 abtol = lsode_get_atol(blsys);
908 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 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
912 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for lsode.");
913 lsodesys.status = lsode_nok;
914 return 4;
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 mf = 21; /* 21 = BDF with exact jacobian. 22 = BDF with finite diff Jacobian */
929
930 if(x[0] > integrator_getsample(blsys, 2)){
931 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value");
932 return 5;
933 }
934
935 /* 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 /* -- store the initial values of all the stuff */
941 integrator_output_write(blsys);
942 integrator_output_write_obs(blsys);
943
944 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
953 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 asc_assert(xend > xprev);
958 /* CONSOLE_DEBUG("LSODE call #%lu: x = [%f,%f]", index,xprev,xend); */
959
960 # ifndef NO_SIGNAL_TRAPS
961 if (SETJMP(g_fpe_env)==0) {
962 # endif /* NO_SIGNAL_TRAPS */
963
964 /* CONSOLE_DEBUG("Calling LSODE with end-time = %f",xend); */
965 /*
966 switch(mf){
967 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 return 0; * failure *
980 }
981 */
982
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 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
993 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 6;
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 lsode_write_istate(istate);
1017 FPRINTF(ASCERR, "\nFurthest point reached was t = %g.\n",x[0]);
1018 error_reporter_end_flush();
1019
1020 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1021 integrator_output_close(blsys);
1022 return 7;
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 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1028 lsodesys.status = lsode_ok; /* clean up before we go */
1029 lsodesys.lastcall = lsode_none;
1030 integrator_output_close(blsys);
1031 return 8;
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 reporterstatus = integrator_output_write(blsys);
1041
1042 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 9;
1049 }
1050
1051 if (nobs > 0) {
1052 # ifndef NO_SIGNAL_TRAPS
1053 if (SETJMP(g_fpe_env)==0) {
1054 # 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 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1073 lsodesys.status = lsode_ok; /* clean up before we go */
1074 lsodesys.lastcall = lsode_none;
1075 integrator_output_close(blsys);
1076 return 10;
1077 }
1078 # endif /* NO_SIGNAL_TRAPS */
1079 }
1080 /* CONSOLE_DEBUG("Integration completed from %3g to %3g.",xprev,x[0]); */
1081 }
1082
1083 CONSOLE_DEBUG("...");
1084 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 CONSOLE_DEBUG("...");
1088
1089
1090 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1091
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 CONSOLE_DEBUG("--- LSODE done ---");
1102 return 0; /* success */
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 11;
1109
1110 #endif
1111 }
1112
1113 /**
1114 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
1118 @see the comments for 'xerrwv' from lsode.f, with which XASCWV is compatible...
1119
1120 @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 0 or 1 means recoverable (control returns to caller).
1125 2 means fatal (run is aborted--see note below).
1126 @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 int *level,
1135 int *ni,
1136 int *i1,
1137 int *i2,
1138 int *nr,
1139 double *r1,
1140 double *r2
1141 ){
1142 static double r1last;
1143
1144 asc_assert(*level!=2); // LSODE doesn't give level 2 in our version.
1145
1146 switch(*nerr){
1147 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 case 204:
1156 if(*nr==0 && *ni==0)return;
1157 if(*nr==2){
1158 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);
1159 return;
1160 } break;
1161 case 205:
1162 if(*nr==2){
1163 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);
1164 return;
1165 } break;
1166 case 27:
1167 if(*nr==1 && *ni==1){
1168 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Trouble with INTDY: itask = %d, tout = %f", *i1, *r1);
1169 return;
1170 } break;
1171 }
1172
1173 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1174
1175 /* note that %.*s means that a string length (integer) and string pointer are being required */
1176 FPRINTF(stderr,"LSODE error: (%d) %.*s",*nerr,*nmes,msg);
1177 if (*ni == 1) {
1178 FPRINTF(stderr,"\nwhere i1 = %d",*i1);
1179 }
1180 if (*ni == 2) {
1181 FPRINTF(stderr,"\nwhere i1 = %d, i2 = %d",*i1,*i2);
1182 }
1183 if (*nr == 1) {
1184 FPRINTF(stderr,"\nwhere r1 = %.13g", *r1);
1185 }
1186 if (*nr == 2) {
1187 FPRINTF(stderr,"\nwhere r1 = %.13g, r2 = %.13g", *r1,*r2);
1188 }
1189 error_reporter_end_flush();
1190 }

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