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

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