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

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