/[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 855 - (show annotations) (download) (as text)
Wed Sep 20 14:00:41 2006 UTC (14 years ago) by johnpye
File MIME type: text/x-csrc
File size: 29952 byte(s)
Removed some verbage in Integrator output
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
59 #include "slv_types.h"
60 #include "mtx.h"
61 #include "rel.h"
62 #include "var.h"
63 #include "discrete.h"
64 #include "conditional.h"
65 #include "bnd.h"
66 #include "logrel.h"
67 #include "slv_common.h"
68 #include "linsol.h"
69 #include "linsolqr.h"
70 #include "slv_client.h"
71
72 #include "integrator.h"
73 #include "lsode.h"
74
75 /*
76 #include "Sensitivity.h"
77 *//* see the packages dir */
78
79 #ifndef lint
80 static CONST char LsodeID[] = "$Id: Lsode.c,v 1.29 2000/01/25 02:26:31 ballan Exp $";
81 #endif
82
83
84 /*
85 * NOUNDERBARS --> FORTRAN compiler naming convention for subroutine
86 * is wierd. WIN32/CRAY is treated as special case
87 */
88 #ifdef APOLLO
89 #define NOUNDERBARS TRUE
90 #endif
91 #ifdef _HPUX_SOURCE
92 #define NOUNDERBARS TRUE
93 #endif
94 /* AIX xlf will not suffix an underbar on a symbol
95 * unless xlf is given the ``-qextname'' option
96 */
97 #ifdef _AIX
98 #define NOUNDERBARS TRUE
99 #endif
100
101 #ifdef NOUNDERBARS
102 #define LSODE lsode
103 #define LSODE_JEX jex
104 #define LSODE_FEX fex
105 #define GETCOMMON get_lsode_common
106 #define XASCWV xascwv
107 #else
108 /* sun, __alpha, __sgi, ... */
109 #define LSODE lsode_
110 #define LSODE_JEX jex_
111 #define LSODE_FEX fex_
112 #define GETCOMMON get_lsode_common_
113 #define XASCWV xascwv_
114 #endif
115
116 #if defined(CRAY) || (defined(__WIN32__) && !defined(__MINGW32_VERSION))
117 #undef LSODE
118 #undef LSODE_JEX
119 #undef LSODE_FEX
120 #undef GETCOMMON
121 #undef XASCWV
122 #define XASCWV XASCWV
123 #define LSODE LSODE
124 #define LSODE_JEX JEX
125 #define LSODE_FEX FEX
126 #define GETCOMMON GET_LSODE_COMMON
127 #endif
128
129 #define DOTIME FALSE
130
131 /* definitions of lsode supported children of atoms, etc */
132 /********************************************************************/
133 /* default input tolerances for lsode */
134 #define RTOLDEF 1e-6
135 #define ATOLDEF 1e-6
136 /* solver_var children expected for state variables */
137 static symchar *g_symbols[2];
138 #define STATERTOL g_symbols[0]
139 #define STATEATOL g_symbols[1]
140 static
141 void InitTolNames(void)
142 {
143 STATERTOL = AddSymbol("ode_rtol");
144 STATEATOL = AddSymbol("ode_atol");
145 }
146
147 /**
148 Because LSODE doesn't seem to make an allowance for 'client data' we
149 have to store this as a 'local global' and fish it out when we're in the
150 callbacks.
151 */
152 IntegratorSystem *l_lsode_blsys;
153
154 enum Lsode_enum {
155 lsode_none, /* true on first call */
156 lsode_function, lsode_derivative, /* functions or gradients done */
157 lsode_sparse, lsode_dense, /* what type of backend should we */
158 lsode_band, /* use for the integrator */
159 lsode_ok, lsode_nok /* bad return from func or grad */
160 };
161
162 static struct Lsode_Data {
163 enum Lsode_enum lastcall; /* type of last call; func or grad */
164 enum Lsode_enum status; /* solve status */
165 int partitioned; /* partioned func evals or not */
166 } lsodesys = {lsode_none, lsode_ok, 1};
167
168
169 /*--------------------------
170 Data space for use by LSODE
171 */
172 typedef struct{
173 long n_eqns; /**< dimension of state vector */
174 int *input_indices; /**< vector of state vars indexes */
175 int *output_indices; /**< vector of derivative var indexes */
176 struct var_variable **y_vars; /**< NULL terminated list of states vars */
177 struct var_variable **ydot_vars; /**< NULL terminated list of derivative vars*/
178 struct rel_relation **rlist; /**< NULL terminated list of relevant rels
179 to be differentiated */
180 double **dydx_dx; /**< change in derivatives wrt states
181 I prefer to call this: d(ydot)/dy */
182 } IntegratorLsodeData;
183
184
185 /*----------------------------
186 Function types that LSODE wants to use
187 */
188
189 /**
190 Type of function used to evaluate derivative system.
191 */
192 typedef void (LsodeEvalFn)( int * ,double * ,double *,double *);
193
194 /**
195 Type of function used to evaluate jacobian system.
196 */
197 typedef void (LsodeJacobianFn)(int *, double *, double *, int *, int *,
198 double *, int *);
199
200 /*----------------------------
201 forward declarations
202 */
203
204 int integrator_lsode_setup_diffs(IntegratorSystem *blsys);
205 static double **MakeDenseMatrix(int nrows, int ncols);
206 static void DestroyDenseMatrix(double **matrix,int nrows);
207
208 /**
209 void LSODE(&fex, &neq, y, &x, &xend, &itol, reltol, abtol, &itask,
210 &istate, &iopt ,rwork, &lrw, iwork, &liw, &jex, &mf);
211
212 This is a prototype for the *fortran* LSODE function.
213
214 No 'extern' here, so we want linker to complain if no static linkage.
215 */
216 void LSODE(LsodeEvalFn*,int *neq ,double *y ,double *x
217 ,double *xend
218 ,int *itol ,double *reltol ,double *abtol
219 ,int *itask ,int *istate ,int *iopt
220 ,double *rwork ,int *lrw
221 ,int *iwork ,int *liw
222 ,LsodeJacobianFn *jex ,int *mf
223 );
224
225 /*------------------------------------------------------
226 Memory allocation/free
227 */
228
229 void integrator_lsode_create(IntegratorSystem *blsys){
230 IntegratorLsodeData *d;
231 d = ASC_NEW_CLEAR(IntegratorLsodeData);
232 d->n_eqns=0;
233 d->input_indices=NULL;
234 d->output_indices=NULL;
235 d->y_vars=NULL;
236 d->ydot_vars=NULL;
237 d->rlist=NULL;
238 d->dydx_dx=NULL;
239 blsys->enginedata=(void*)d;
240 }
241
242 /**
243 Cleanup the data struct that belongs to BLSODE
244 */
245 void integrator_lsode_free(void *enginedata){
246 IntegratorLsodeData d;
247 d = *((IntegratorLsodeData *)enginedata);
248
249 if(d.input_indices)ASC_FREE(d.input_indices);
250 d.input_indices = NULL;
251
252 if(d.output_indices)ASC_FREE(d.output_indices);
253 d.output_indices = NULL;
254
255 if(d.y_vars)ASC_FREE(d.y_vars);
256 d.y_vars = NULL;
257
258 if(d.ydot_vars)ASC_FREE(d.ydot_vars);
259 d.ydot_vars = NULL;
260
261 if(d.rlist)ASC_FREE(d.rlist);
262 d.rlist = NULL;
263
264 if(d.dydx_dx)DestroyDenseMatrix(d.dydx_dx, d.n_eqns);
265 d.dydx_dx = NULL;
266
267 d.n_eqns = 0L;
268 }
269
270 /*-------------------------------------------------------
271 Prototype for the LSODE function (c.f. lsode.f)
272 */
273
274 /*---------------------------------------------------------
275 Couple of matrix methods...?
276 */
277
278 static double **MakeDenseMatrix(int nrows, int ncols){
279 int c;
280 double **result;
281 assert(nrows>0);
282 assert(ncols>0);
283 result = ASC_NEW_ARRAY(double *, nrows);
284 for (c=0;c<nrows;c++) {
285 result[c] = ASC_NEW_ARRAY_CLEAR(double, ncols);
286 }
287 return result;
288 }
289
290 static void DestroyDenseMatrix(double **matrix,int nrows){
291 int c;
292 if (matrix) {
293 for (c=0;c<nrows;c++) {
294 if (matrix[c]) {
295 ascfree((char *)matrix[c]);
296 }
297 }
298 ascfree((char *)matrix);
299 }
300 }
301
302 /*------------------------------------------------------------------------------
303 PROBLEM ANALYSIS
304 */
305
306 /**
307 @TODO needs work. Assumes struct Instance* and struct var_variable*
308 are synonymous, which demonstrates the need for a method to take
309 an instance and ask the solvers for its global or local index
310 if var and inst are decoupled.
311 */
312 int integrator_lsode_setup_diffs(IntegratorSystem *blsys) {
313 long n_eqns;
314 unsigned long nch,i;
315
316 struct var_variable **vp;
317 int *ip;
318
319 IntegratorLsodeData *enginedata;
320 enginedata = (IntegratorLsodeData *)blsys->enginedata;
321 assert(enginedata!=NULL);
322
323 assert(enginedata->n_eqns==blsys->n_y);
324
325 /*
326 Put the
327 Let us now process what we consider *inputs* to the problem as
328 far as ASCEND is concerned; i.e. the state vars or the y_vars's
329 if you prefer.
330 */
331 nch = enginedata->n_eqns;
332
333
334 vp = enginedata->y_vars;
335 ip = enginedata->input_indices;
336 for (i=0;i<nch;i++) {
337 *vp = (struct var_variable *)blsys->y[i];
338 *ip = var_sindex(*vp);
339 vp++;
340 ip++;
341 }
342 *vp = NULL; /* terminate */
343
344 /*
345 Let us now go for the outputs, ie the derivative terms.
346 */
347 vp = enginedata->ydot_vars;
348 ip = enginedata->output_indices;
349 for (i=0;i<nch;i++) {
350 *vp = (struct var_variable *)blsys->ydot[i];
351 *ip = var_sindex(*vp);
352 vp++; /* dont assume that a var is synonymous with */
353 ip++; /* an Instance; that might/will change soon */
354 }
355 *vp = NULL; /* terminate */
356
357 return 0;
358 }
359
360
361
362 /**
363 allocates, fills, and returns the atol vector based on BLSODE
364
365 State variables missing child ode_rtol will be defaulted to ATOLDEF
366 */
367 static double *blsode_get_atol( IntegratorSystem *blsys) {
368
369 struct Instance *tol;
370 double *atoli;
371 int i,len;
372
373 len = blsys->n_y;
374 atoli = ASC_NEW_ARRAY(double, blsys->n_y+1);
375 if (atoli == NULL) {
376 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
377 return atoli;
378 }
379 InitTolNames();
380 for (i=0; i<len; i++) {
381 tol = ChildByChar(var_instance(blsys->y[i]),STATEATOL);
382 if (tol == NULL || !AtomAssigned(tol) ) {
383 atoli[i] = ATOLDEF;
384 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming atol = %3g"
385 "for ode_atol child undefined for state variable %ld."
386 ,ATOLDEF, blsys->y_id[i]
387 );
388 } else {
389 atoli[i] = RealAtomValue(tol);
390 }
391 }
392 atoli[len] = ATOLDEF;
393 return atoli;
394 }
395
396 /**
397 Allocates, fills, and returns the rtol vector based on BLSODE
398
399 State variables missing child ode_rtol will be defaulted to RTOLDEF
400 */
401 static double *blsode_get_rtol( IntegratorSystem *blsys) {
402
403 struct Instance *tol;
404 double *rtoli;
405 int i,len;
406
407 len = blsys->n_y;
408 rtoli = ASC_NEW_ARRAY(double, blsys->n_y+1);
409 if (rtoli == NULL) {
410 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
411 return rtoli;
412 }
413 InitTolNames();
414 for (i=0; i<len; i++) {
415 tol = ChildByChar(var_instance(blsys->y[i]),STATERTOL);
416 if (tol == NULL || !AtomAssigned(tol) ) {
417 rtoli[i] = RTOLDEF;
418
419 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming rtol = %3g"
420 "for ode_rtol child undefined for state variable %ld."
421 ,ATOLDEF, blsys->y_id[i]
422 );
423
424 } else {
425 rtoli[i] = RealAtomValue(tol);
426 }
427 }
428 rtoli[len] = RTOLDEF;
429 return rtoli;
430 }
431
432 /*
433 Write out a a status message based on the istate parameter.
434 */
435 static void blsode_write_istate( int istate) {
436 switch (istate) {
437 case -1:
438 FPRINTF(ASCERR,"Excess steps taken on this call (perhaps wrong MF).");
439 break;
440 case -2:
441 FPRINTF(ASCERR,"Excess accuracy requested (tolerances too small).");
442 break;
443 case -3:
444 FPRINTF(ASCERR,"Illegal input detected.");
445 break;
446 case -4:
447 FPRINTF(ASCERR,"Repeated error test failures.");
448 break;
449 case -5:
450 FPRINTF(ASCERR,"Repeated convergence failures.");
451 break;
452 case -6:
453 FPRINTF(ASCERR,"Error weight became zero during problem.");
454 break;
455 case -7:
456 FPRINTF(ASCERR,"User patience became zero during problem.");
457 break;
458 default:
459 FPRINTF(ASCERR,"Unknown error code %d.",istate);
460 break;
461 }
462 }
463
464 /**
465 Free memory allocated for the LSODE, but first check.
466 */
467 static void blsode_free_mem(double *y, double *reltol, double *abtol, double *rwork,
468 int *iwork, double *obs, double *dydx)
469 {
470 if (y != NULL) {
471 ascfree((double *)y);
472 }
473 if (reltol != NULL) {
474 ascfree((double *)reltol);
475 }
476 if (abtol != NULL) {
477 ascfree((double *)abtol);
478 }
479 if (rwork != NULL) {
480 ascfree((double *)rwork);
481 }
482 if (iwork != NULL) {
483 ascfree((int *)iwork);
484 }
485 if (obs != NULL) {
486 ascfree((double *)obs);
487 }
488 if (dydx != NULL) {
489 ascfree((double *)dydx);
490 }
491 }
492
493 /*
494 *********************************************************************
495 * This code is provided for the benefit of a temporary
496 * fix for the derivative problem in Lsode.
497 * The proper permanent fix for lsode is to dump it in favor of
498 * cvode or dassl.
499 * Extended 7/95 baa to deal with linsolqr and blsode.
500 * It is assumed the system has been solved at the current point.
501 *********************************************************************
502 */
503 int Asc_BLsodeDerivatives(slv_system_t sys, double **dy_dx,
504 int *inputs_ndx_list, int ninputs,
505 int *outputs_ndx_list, int noutputs)
506 {
507 static int n_calls = 0;
508 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
509 mtx_matrix_t mtx;
510 int32 capacity;
511 real64 *scratch_vector = NULL;
512 int result=0;
513
514 (void)NumberFreeVars(NULL); /* used to re-init the system */
515 (void)NumberIncludedRels(NULL); /* used to re-init the system */
516 if (!sys) {
517 FPRINTF(stderr,"The solve system does not exist !\n");
518 return 1;
519 }
520
521 result = Compute_J(sys);
522 if (result) {
523 FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
524 return 1;
525 }
526
527 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
528 if (lqr_sys==NULL) {
529 FPRINTF(stderr,"Early termination due to missing linsolqr system.\n");
530 return 1;
531 }
532 mtx = slv_get_sys_mtx(sys); /* get the matrix */
533 if (mtx==NULL) {
534 FPRINTF(stderr,"Early termination due to missing mtx in linsolqr.\n");
535 return 1;
536 }
537 capacity = mtx_capacity(mtx);
538 scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
539 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
540
541 result = LUFactorJacobian(sys);
542 if (result) {
543 FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
544 goto error;
545 }
546 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
547 inputs_ndx_list, ninputs,
548 outputs_ndx_list, noutputs);
549
550 linsolqr_remove_rhs(lqr_sys,scratch_vector);
551 if (result) {
552 FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
553 goto error;
554 }
555
556 error:
557 n_calls++;
558 if (scratch_vector) {
559 ascfree((char *)scratch_vector);
560 }
561 return result;
562 }
563
564 /**
565 The current way that we are getting the derivatives (if the problem
566 was solved partitioned) messes up the slv_system so that we *have*
567 to do a *presolve* rather than a simply a *resolve* before doing
568 function calls. This code below attempts to handle these cases.
569 */
570 static void LSODE_FEX( int *n_eq ,double *t ,double *y ,double *ydot)
571 {
572 slv_status_t status;
573
574 /* slv_parameters_t parameters; pity lsode doesn't allow error returns */
575 int i;
576 unsigned long ok;
577
578 #if DOTIME
579 double time1,time2;
580 #endif
581
582 #if DOTIME
583 CONSOLE_DEBUG("Calling for a function evaluation");
584 time1 = tm_cpu_time();
585 #endif
586
587 /*
588 t[1]=t[0]; can't do this. lsode calls us with a different t than the x we sent in.
589 */
590 integrator_set_t(l_lsode_blsys, t[0]);
591 integrator_set_y(l_lsode_blsys, y);
592
593 #if DOTIME
594 time2 = tm_cpu_time();
595 #endif
596
597 switch(lsodesys.lastcall) {
598 case lsode_none: /* first call */
599 case lsode_derivative:
600 if (lsodesys.partitioned) {
601 slv_presolve(l_lsode_blsys->system);
602 } else {
603 slv_resolve(l_lsode_blsys->system);
604 }
605 break;
606 default:
607 case lsode_function:
608 slv_resolve(l_lsode_blsys->system);
609 break;
610 }
611 slv_solve(l_lsode_blsys->system);
612 slv_get_status(l_lsode_blsys->system, &status);
613 ok = integrator_checkstatus(status);
614
615 #if DOTIME
616 time2 = tm_cpu_time() - time2;
617 #endif
618
619 if (!ok) {
620 ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
621 FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");
622 for (i = 0; i< *n_eq; i++) {
623 FPRINTF(ASCERR,"y[%4d] = %f\n",i, y[i]);
624 }
625 error_reporter_end_flush();
626 lsodesys.status = lsode_nok;
627 } else {
628 lsodesys.status = lsode_ok;
629 }
630 integrator_get_ydot(l_lsode_blsys, ydot);
631
632 lsodesys.lastcall = lsode_function;
633 #if DOTIME
634 time1 = tm_cpu_time() - time1;
635 CONSOLE_DEBUG("Function evalulation has been completed in time %g. True function call time = %g",time1,time2);
636 #endif
637 }
638
639 /**
640 Evaluate the jacobian
641 */
642 static void LSODE_JEX(int *neq ,double *t, double *y,
643 int *ml ,int *mu ,double *pd, int *nrpd)
644 {
645 int nok = 0;
646 int i,j;
647
648 IntegratorLsodeData enginedata=*((IntegratorLsodeData *)l_lsode_blsys->enginedata);
649
650 UNUSED_PARAMETER(t);
651 UNUSED_PARAMETER(y);
652 UNUSED_PARAMETER(ml);
653 UNUSED_PARAMETER(mu);
654
655 #if DOTIME
656 double time1;
657
658 CONSOLE_DEBUG("Calling for a gradient evaluation");
659 time1 = tm_cpu_time();
660 #endif
661 /*
662 * Make the real call.
663 */
664 nok = Asc_BLsodeDerivatives(l_lsode_blsys->system
665 , enginedata.dydx_dx
666 , enginedata.input_indices
667 , *neq
668 , enginedata.output_indices
669 , *nrpd
670 );
671
672 if (nok) {
673 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in computing the derivatives for the system. Failing...");
674 lsodesys.status = lsode_nok;
675 lsodesys.lastcall = lsode_derivative;
676 return;
677 } else {
678 lsodesys.status = lsode_ok;
679 lsodesys.lastcall = lsode_derivative;
680 }
681 /*
682 * Map data from C based matrix to Fortan matrix.
683 * We will send in a column major ordering vector for pd.
684 */
685 for (j=0;j<*neq;j++) {
686 for (i=0;i<*nrpd;i++) { /* column wise - hence rows change faster */
687 *pd++ = enginedata.dydx_dx[i][j];
688 }
689 }
690
691 #if DOTIME
692 time1 = tm_cpu_time() - time1;
693 CONSOLE_DEBUG("Time to do gradient evaluation %g",time1);
694 #endif
695
696 return;
697 }
698
699 /**
700 The public function: here we do the actual integration, I guess.
701 */
702 int integrator_lsode_solve(IntegratorSystem *blsys
703 , unsigned long start_index, unsigned long finish_index
704 ){
705 slv_status_t status;
706 slv_parameters_t params;
707 IntegratorLsodeData *d;
708
709 double x[2];
710 double xend,xprev;
711 unsigned long nsamples, neq;
712 long nobs;
713 int itol, itask, mf, lrw, liw;
714 unsigned long index;
715 int istate, iopt;
716 double * rwork;
717 int * iwork;
718 double *y, *abtol, *reltol, *obs, *dydx;
719 int my_neq;
720 FILE *y_out =NULL;
721 FILE *obs_out =NULL;
722
723 /* store the local variable so that we can get at stuff from inside LSODE_FEX. */
724 l_lsode_blsys = blsys;
725
726 d = (IntegratorLsodeData *)(blsys->enginedata);
727
728 /* the numer of equations must be equal to blsys->n_y, the number of states */
729 d->n_eqns = blsys->n_y;
730 assert(d->n_eqns>0);
731
732 d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
733 d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
734 d->dydx_dx = MakeDenseMatrix(d->n_eqns,d->n_eqns);
735
736 d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);
737 d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);
738
739 integrator_lsode_setup_diffs(blsys);
740
741 /* this is a lie, but we will keep it.
742 We handle any linsol/linsolqr based solver. */
743 if (strcmp(slv_solver_name(slv_get_selected_solver(blsys->system)),"QRSlv") != 0) {
744 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv must be selected before integration.");
745 return 0;
746 }
747
748 slv_get_status(l_lsode_blsys->system, &status);
749
750 if (status.struct_singular) {
751 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system is structurally singular.");
752 lsodesys.status = lsode_nok;
753 return 0;
754 }
755
756 #if defined(STATIC_LSOD) || defined (DYNAMIC_LSOD)
757
758 /* here we assume integrators.c is in charge of dynamic loading */
759
760 slv_get_parameters(blsys->system,&params);
761 lsodesys.partitioned = 1;
762
763 nsamples = integrator_getnsamples(blsys);
764 if (nsamples <2) {
765 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined.");
766 lsodesys.status = lsode_nok;
767 return 0;
768 }
769 neq = blsys->n_y;
770 nobs = blsys->n_obs;
771
772 x[0] = integrator_get_t(blsys);
773 x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */
774 lrw = 22 + 9*neq + neq*neq;
775 rwork = ASC_NEW_ARRAY_CLEAR(double, lrw+1);
776 liw = 20 + neq;
777 iwork = ASC_NEW_ARRAY_CLEAR(int, liw+1);
778 y = integrator_get_y(blsys, NULL);
779 reltol = blsode_get_rtol(blsys);
780 abtol = blsode_get_atol(blsys);
781 obs = integrator_get_observations(blsys, NULL);
782 dydx = ASC_NEW_ARRAY_CLEAR(double, neq+1);
783 if (!y || !obs || !abtol || !reltol || !rwork || !iwork || !dydx) {
784 blsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
785 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for blsode.");
786 lsodesys.status = lsode_nok;
787 return 0;
788 }
789
790 /*
791 Prepare args and call lsode.
792 */
793 itol = 4;
794 itask = 1;
795 istate = 1;
796 iopt = 1;
797 rwork[4] = integrator_get_stepzero(blsys);
798 rwork[5] = integrator_get_maxstep(blsys);
799 rwork[6] = integrator_get_minstep(blsys);
800 iwork[5] = integrator_get_maxsubsteps(blsys);
801 mf = 21; /* gradients 22 -- implies finite diff */
802
803 /* put the values from derivative system into the record */
804 integrator_setsample(blsys, start_index, x[0]);
805
806 integrator_output_init(blsys);
807
808 my_neq = (int)neq;
809
810 /*
811 First time entering lsode, x is input. After that,
812 lsode uses x as output (y output is y(x)). To drive
813 the loop ahead in time, all we need to do is keep upping
814 xend.
815 */
816 blsys->currentstep = 0;
817 for (index = start_index; index < finish_index; index++, blsys->currentstep++) {
818 xend = integrator_getsample(blsys, index+1);
819 xprev = x[0];
820 /* CONSOLE_DEBUG("BEFORE %lu BLSODE CALL\n", index); */
821
822 # ifndef NO_SIGNAL_TRAPS
823 if (setjmp(g_fpe_env)==0) {
824 # endif /* NO_SIGNAL_TRAPS */
825
826
827 LSODE(&(LSODE_FEX), &my_neq, y, x, &xend,
828 &itol, reltol, abtol, &itask, &istate,
829 &iopt ,rwork, &lrw, iwork, &liw, &(LSODE_JEX), &mf);
830
831
832 # ifndef NO_SIGNAL_TRAPS
833 } else {
834 FPRINTF(stderr,
835 "Integration terminated due to float error in LSODE call.\n");
836 blsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
837 lsodesys.status = lsode_ok; /* clean up before we go */
838 lsodesys.lastcall = lsode_none;
839 if (y_out!=NULL) {
840 fclose(y_out);
841 }
842 if (obs_out!=NULL) {
843 fclose(obs_out);
844 }
845 return 0;
846 }
847 # endif /* NO_SIGNAL_TRAPS */
848
849 /* CONSOLE_DEBUG("AFTER %lu LSODE CALL\n", index); */
850 /* this check is better done in fex,jex, but lsode takes no status */
851 if (Solv_C_CheckHalt()) {
852 if (istate >= 0) {
853 istate=-7;
854 }
855 }
856
857 if (istate < 0 ) {
858 /* some kind of error occurred... */
859 ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
860 FPRINTF(ASCERR, "In LSODE integrator: index = %lu, ", index);
861 FPRINTF(ASCERR, "istate = %d.\n", istate);
862 FPRINTF(ASCERR, "Farthest point reached: t = %g.\n",x[0]);
863 blsode_write_istate(istate);
864 error_reporter_end_flush();
865
866 blsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
867 integrator_output_close(blsys);
868 return 0;
869 }
870
871 if (lsodesys.status==lsode_nok) {
872 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to an error in derivative computations.");
873 blsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
874 lsodesys.status = lsode_ok; /* clean up before we go */
875 lsodesys.lastcall = lsode_none;
876 integrator_output_close(blsys);
877 return 0;
878 }
879
880 integrator_setsample(blsys, index+1, x[0]);
881 /* record when lsode actually came back */
882 integrator_set_t(blsys, x[0]);
883 integrator_set_y(blsys, y);
884 /* put x,y in d in case lsode got x,y by interpolation, as it does */
885
886 integrator_output_write(blsys);
887
888 if (nobs > 0) {
889 # ifndef NO_SIGNAL_TRAPS
890 if (setjmp(g_fpe_env)==0) {
891 # endif /* NO_SIGNAL_TRAPS */
892
893 /* solve for obs since d isn't necessarily already
894 computed there though lsode's x and y may be.
895 Note that since lsode usually steps beyond xend
896 x1 usually wouldn't be x0 precisely if the x1/x0
897 scheme worked, which it doesn't anyway. */
898
899 LSODE_FEX(&my_neq, x, y, dydx);
900
901 /* calculate observations, if any, at returned x and y. */
902 obs = integrator_get_observations(blsys, obs);
903
904 integrator_output_write_obs(blsys);
905
906 # ifndef NO_SIGNAL_TRAPS
907 } else {
908 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in LSODE FEX call.");
909 blsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
910 lsodesys.status = lsode_ok; /* clean up before we go */
911 lsodesys.lastcall = lsode_none;
912 integrator_output_close(blsys);
913 return 0;
914 }
915 # endif /* NO_SIGNAL_TRAPS */
916 }
917 /* CONSOLE_DEBUG("Integration completed from %3g to %3g.",xprev,x[0]); */
918 }
919
920 CONSOLE_DEBUG("");
921 CONSOLE_DEBUG("Number of steps taken: %1d.", iwork[10]);
922 CONSOLE_DEBUG("Number of function evaluations: %1d.", iwork[11]);
923 CONSOLE_DEBUG("Number of Jacobian evaluations: %1d.", iwork[12]);
924 CONSOLE_DEBUG("");
925
926
927 blsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
928
929 /*
930 * return the system to its original state.
931 */
932
933 lsodesys.status = lsode_ok;
934 lsodesys.lastcall = lsode_none;
935
936 integrator_output_close(blsys);
937
938 CONSOLE_DEBUG("--- BLSODE done ---");
939 return 1;
940
941 #else /* STATIC_LSOD || DYNAMIC_LSOD */
942
943 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration will not be performed. LSODE binary not available.");
944 lsodesys.status = lsode_nok;
945 return 0;
946
947 #endif
948 }
949
950 /**
951 This function is XASCWV, a compatible
952 replacement for the function XERRWV from lsode.f by Alan Hindmarsh.
953
954 To the extent permitted by the GPL that covers ASCEND, the following
955 function is placed in the public domain, because LSODE itself was
956 released in the public domain.
957
958 The following is the FORTRAN header for 'xerrwv' from LSODE,
959 with which XASCWV is compatible...
960
961 *-----------------------------------------------------------------------
962 Subroutine XERRWV
963 A simplified version of the slatec error handling package.
964 Written by A. C. Hindmarsh at LLNL. Version of March 30, 1987.
965 This version is in double precision.
966
967 All arguments are input arguments.
968
969 msg = the message (hollerith literal or integer array).
970 nmes = the length of msg (number of characters).
971 nerr = the error number (not used).
972 level = the error level..
973 0 or 1 means recoverable (control returns to caller).
974 2 means fatal (run is aborted--see note below).
975 ni = number of integers (0, 1, or 2) to be printed with message.
976 i1,i2 = integers to be printed, depending on ni.
977 nr = number of reals (0, 1, or 2) to be printed with message.
978 r1,r2 = reals to be printed, depending on nr.
979
980 NOTE.. This routine is machine-dependent and specialized for use
981 in limited context, in the following ways..
982 1. The number of hollerith characters stored per word, denoted
983 by ncpw below, is a data-loaded constant.
984 2. The value of nmes is assumed to be at most 60.
985 (multi-line messages are generated by repeated calls.)
986 3. If level = 2, control passes to the statement stop
987 to abort the run. this statement may be machine-dependent.
988 4. r1 and r2 are assumed to be in double precision and are printed
989 in d21.13 format.
990 5. The common block /eh0001/ below is data-loaded (a machine-
991 dependent feature) with default values.
992 This block is needed for proper retention of parameters used by
993 this routine which the user can reset by calling xsetf or xsetun.
994 The variables in this block are as follows..
995
996 mesflg = print control flag..
997 1 means print all messages (the default).
998 0 means no printing.
999 lunit = logical unit number for messages.
1000 the default is 6 (machine-dependent).
1001 *-----------------------------------------------------------------------
1002 The following are instructions for installing this routine
1003 in different machine environments.
1004
1005 To change the default output unit, change the data statement
1006 in the block data subprogram below.
1007
1008 For a different number of characters per word, change the
1009 data statement setting ncpw below, and format 10. alternatives for
1010 various computers are shown in comment cards.
1011
1012 For a different run-abort command, change the statement following
1013 statement 100 at the end.
1014 *-----------------------------------------------------------------------
1015 */
1016 void XASCWV( char *msg, /* array of char/int not NULL ended, len *nmes */
1017 int *nmes, /* len of msg */
1018 int *nerr,
1019 int *level,
1020 int *ni,
1021 int *i1,
1022 int *i2,
1023 int *nr,
1024 double *r1,
1025 double *r2
1026 )
1027 {
1028 (void)nerr;
1029
1030 /* write the message. lot easier in C, geez */
1031 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1032
1033 /* note that %.*s means that a string length (integer) and string pointer are being required */
1034 FPRINTF(stderr,"LSODE error: %.*s",*nmes,msg);
1035 if (*ni == 1) {
1036 FPRINTF(stderr,"\nwhere i1 = %d\n",*i1);
1037 }
1038 if (*ni == 2) {
1039 FPRINTF(stderr,"\nwhere i1 = %d, i2 = %d",*i1,*i2);
1040 }
1041 if (*nr == 1) {
1042 FPRINTF(stderr,"\nwhere r1 = %.13g", *r1);
1043 }
1044 if (*nr == 2) {
1045 FPRINTF(stderr,"\nwhere r1 = %.13g, r2 = %.13g", *r1,*r2);
1046 }
1047 if (*level != 2) {
1048 error_reporter_end_flush();
1049 return;
1050 }
1051
1052 /* NOT reached. lsode does NOT make level 2 calls in our version. */
1053 error_reporter_end_flush();
1054 Asc_Panic(3,"xascwv", "LSODE really really confused");
1055
1056 return; /* not reached */
1057 }

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