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

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