/[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 942 - (show annotations) (download) (as text)
Sat Nov 25 05:26:47 2006 UTC (14 years, 2 months ago) by johnpye
File MIME type: text/x-csrc
File size: 30912 byte(s)
Incorporated 'SolverParameters' functionality into 'Integrator', both at C level and C++/Python.
This shouldn't break Tcl/Tk as default parameters will be set and used invisibly.
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 PARAMETERS
272 */
273
274 enum ida_parameters{
275 LSODE_PARAM_TIMING
276 ,LSODE_PARAMS_SIZE
277 };
278
279 /**
280 Here the full set of parameters is defined, along with upper/lower bounds,
281 etc. The values are stuck into the blsys->params structure.
282
283 @return 0 on success
284 */
285 int integrator_lsode_params_default(IntegratorSystem *blsys){
286
287 asc_assert(blsys!=NULL);
288 asc_assert(blsys->engine==INTEG_IDA);
289 slv_parameters_t *p;
290 p = &(blsys->params);
291
292 slv_destroy_parms(p);
293
294 if(p->parms==NULL){
295 CONSOLE_DEBUG("params NULL");
296 p->parms = ASC_NEW_ARRAY(struct slv_parameter, LSODE_PARAMS_SIZE);
297 if(p->parms==NULL)return -1;
298 p->dynamic_parms = 1;
299 }else{
300 asc_assert(p->num_parms == LSODE_PARAMS_SIZE);
301 CONSOLE_DEBUG("reusing parm memory");
302 }
303
304 /* reset the number of parameters to zero so that we can check it at the end */
305 p->num_parms = 0;
306
307 slv_param_bool(p,LSODE_PARAM_TIMING
308 ,(SlvParameterInitBool){{"timing"
309 ,"Output timing statistics?",1,NULL
310 }, TRUE}
311 );
312
313 asc_assert(p->num_parms == LSODE_PARAMS_SIZE);
314
315 CONSOLE_DEBUG("Created %d params", p->num_parms);
316
317 return 0;
318 }
319
320 /*---------------------------------------------------------
321 Couple of matrix methods...?
322 */
323
324 static double **lsode_densematrix_create(int nrows, int ncols){
325 int c;
326 double **result;
327 assert(nrows>0);
328 assert(ncols>0);
329 result = ASC_NEW_ARRAY(double *, nrows);
330 for (c=0;c<nrows;c++) {
331 result[c] = ASC_NEW_ARRAY_CLEAR(double, ncols);
332 }
333 return result;
334 }
335
336 static void lsode_densematrix_destroy(double **matrix,int nrows){
337 int c;
338 if (matrix) {
339 for (c=0;c<nrows;c++) {
340 if (matrix[c]) {
341 ascfree((char *)matrix[c]);
342 }
343 }
344 ascfree((char *)matrix);
345 }
346 }
347
348 /*------------------------------------------------------------------------------
349 PROBLEM ANALYSIS
350 */
351
352 /**
353 @TODO needs work. Assumes struct Instance* and struct var_variable*
354 are synonymous, which demonstrates the need for a method to take
355 an instance and ask the solvers for its global or local index
356 if var and inst are decoupled.
357 */
358 int integrator_lsode_setup_diffs(IntegratorSystem *blsys) {
359 /* long n_eqns; */
360 unsigned long nch,i;
361
362 struct var_variable **vp;
363 int *ip;
364
365 IntegratorLsodeData *enginedata;
366 enginedata = (IntegratorLsodeData *)blsys->enginedata;
367 assert(enginedata!=NULL);
368
369 assert(enginedata->n_eqns==blsys->n_y);
370
371 /*
372 Put the
373 Let us now process what we consider *inputs* to the problem as
374 far as ASCEND is concerned; i.e. the state vars or the y_vars's
375 if you prefer.
376 */
377 nch = enginedata->n_eqns;
378
379
380 vp = enginedata->y_vars;
381 ip = enginedata->input_indices;
382 for (i=0;i<nch;i++) {
383 *vp = (struct var_variable *)blsys->y[i];
384 *ip = var_sindex(*vp);
385 vp++;
386 ip++;
387 }
388 *vp = NULL; /* terminate */
389
390 /*
391 Let us now go for the outputs, ie the derivative terms.
392 */
393 vp = enginedata->ydot_vars;
394 ip = enginedata->output_indices;
395 for (i=0;i<nch;i++) {
396 *vp = (struct var_variable *)blsys->ydot[i];
397 *ip = var_sindex(*vp);
398 vp++; /* dont assume that a var is synonymous with */
399 ip++; /* an Instance; that might/will change soon */
400 }
401 *vp = NULL; /* terminate */
402
403 return 0;
404 }
405
406 /**
407 allocates, fills, and returns the atol vector based on LSODE
408
409 State variables missing child ode_rtol will be defaulted to ATOLDEF
410 */
411 static double *lsode_get_atol( IntegratorSystem *blsys) {
412
413 struct Instance *tol;
414 double *atoli;
415 int i,len;
416
417 len = blsys->n_y;
418 atoli = ASC_NEW_ARRAY(double, blsys->n_y+1);
419 if (atoli == NULL) {
420 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
421 return atoli;
422 }
423 InitTolNames();
424 for (i=0; i<len; i++) {
425 tol = ChildByChar(var_instance(blsys->y[i]),STATEATOL);
426 if (tol == NULL || !AtomAssigned(tol) ) {
427 atoli[i] = ATOLDEF;
428 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming atol = %3g"
429 "for ode_atol child undefined for state variable %ld."
430 ,ATOLDEF, blsys->y_id[i]
431 );
432 } else {
433 atoli[i] = RealAtomValue(tol);
434 CONSOLE_DEBUG("Using tolerance %3g for state variable %ld.",atoli[i], blsys->y_id[i]);
435 }
436 }
437 atoli[len] = ATOLDEF;
438 return atoli;
439 }
440
441 /**
442 Allocates, fills, and returns the rtol vector based on LSODE
443
444 State variables missing child ode_rtol will be defaulted to RTOLDEF
445 */
446 static double *lsode_get_rtol( IntegratorSystem *blsys) {
447
448 struct Instance *tol;
449 double *rtoli;
450 int i,len;
451
452 len = blsys->n_y;
453 rtoli = ASC_NEW_ARRAY(double, blsys->n_y+1);
454 if (rtoli == NULL) {
455 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
456 return rtoli;
457 }
458 InitTolNames();
459 for (i=0; i<len; i++) {
460 tol = ChildByChar(var_instance(blsys->y[i]),STATERTOL);
461 if (tol == NULL || !AtomAssigned(tol) ) {
462 rtoli[i] = RTOLDEF;
463
464 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming rtol = %3g"
465 "for ode_rtol child undefined for state variable %ld."
466 ,ATOLDEF, blsys->y_id[i]
467 );
468
469 } else {
470 rtoli[i] = RealAtomValue(tol);
471 }
472 }
473 rtoli[len] = RTOLDEF;
474 return rtoli;
475 }
476
477 /*
478 Write out a a status message based on the istate parameter.
479 */
480 static void lsode_write_istate( int istate) {
481 switch (istate) {
482 case -1:
483 FPRINTF(ASCERR,"Excess steps taken on this call (perhaps wrong MF).");
484 break;
485 case -2:
486 FPRINTF(ASCERR,"Excess accuracy requested (tolerances too small).");
487 break;
488 case -3:
489 FPRINTF(ASCERR,"Illegal input detected (see console).");
490 break;
491 case -4:
492 FPRINTF(ASCERR,"Repeated error test failures (check all inputs).");
493 break;
494 case -5:
495 FPRINTF(ASCERR,"Repeated convergence failures (perhaps bad Jacobian supplied, or wrong choice of MF or tolerances).");
496 break;
497 case -6:
498 FPRINTF(ASCERR,"Error weight became zero during problem (solution component i vanished, and atol or atol(i) = 0).");
499 break;
500 case -7:
501 FPRINTF(ASCERR,"Interrupted? User cancelled operation?");
502 break;
503 default:
504 FPRINTF(ASCERR,"Unknown 'istate' error code %d from LSODE.",istate);
505 break;
506 }
507 }
508
509 /**
510 Free memory allocated for the LSODE, but first check.
511 */
512 static void lsode_free_mem(double *y, double *reltol, double *abtol, double *rwork,
513 int *iwork, double *obs, double *dydx)
514 {
515 if (y != NULL) {
516 ascfree((double *)y);
517 }
518 if (reltol != NULL) {
519 ascfree((double *)reltol);
520 }
521 if (abtol != NULL) {
522 ascfree((double *)abtol);
523 }
524 if (rwork != NULL) {
525 ascfree((double *)rwork);
526 }
527 if (iwork != NULL) {
528 ascfree((int *)iwork);
529 }
530 if (obs != NULL) {
531 ascfree((double *)obs);
532 }
533 if (dydx != NULL) {
534 ascfree((double *)dydx);
535 }
536 }
537
538 /*
539 *********************************************************************
540 * This code is provided for the benefit of a temporary
541 * fix for the derivative problem in Lsode.
542 * The proper permanent fix for lsode is to dump it in favor of
543 * cvode or dassl.
544 * Extended 7/95 baa to deal with linsolqr and lsode.
545 * It is assumed the system has been solved at the current point.
546 *********************************************************************
547 */
548 int lsode_derivatives(slv_system_t sys, double **dy_dx,
549 int *inputs_ndx_list, int ninputs,
550 int *outputs_ndx_list, int noutputs)
551 {
552 static int n_calls = 0;
553 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
554 mtx_matrix_t mtx;
555 int32 capacity;
556 real64 *scratch_vector = NULL;
557 int result=0;
558
559 (void)NumberFreeVars(NULL); /* used to re-init the system */
560 (void)NumberIncludedRels(NULL); /* used to re-init the system */
561 if (!sys) {
562 FPRINTF(stderr,"The solve system does not exist !\n");
563 return 1;
564 }
565
566 result = Compute_J(sys);
567 if (result) {
568 FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
569 return 1;
570 }
571
572 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
573 if (lqr_sys==NULL) {
574 FPRINTF(stderr,"Early termination due to missing linsolqr system.\n");
575 return 1;
576 }
577 mtx = slv_get_sys_mtx(sys); /* get the matrix */
578 if (mtx==NULL) {
579 FPRINTF(stderr,"Early termination due to missing mtx in linsolqr.\n");
580 return 1;
581 }
582 capacity = mtx_capacity(mtx);
583 scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
584 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
585
586 result = LUFactorJacobian(sys);
587 if (result) {
588 FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
589 goto error;
590 }
591 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
592 inputs_ndx_list, ninputs,
593 outputs_ndx_list, noutputs);
594
595 linsolqr_remove_rhs(lqr_sys,scratch_vector);
596 if (result) {
597 FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
598 goto error;
599 }
600
601 error:
602 n_calls++;
603 if (scratch_vector) {
604 ascfree((char *)scratch_vector);
605 }
606 return result;
607 }
608
609 /**
610 The current way that we are getting the derivatives (if the problem
611 was solved partitioned) messes up the slv_system so that we *have*
612 to do a *presolve* rather than a simply a *resolve* before doing
613 function calls. This code below attempts to handle these cases.
614 */
615 static void LSODE_FEX( int *n_eq ,double *t ,double *y ,double *ydot)
616 {
617 slv_status_t status;
618
619 /* slv_parameters_t parameters; pity lsode doesn't allow error returns */
620 /* int i; */
621 unsigned long ok;
622
623 #if DOTIME
624 double time1,time2;
625 #endif
626
627 /* CONSOLE_DEBUG("Calling for a function evaluation"); */
628
629 #if DOTIME
630 CONSOLE_DEBUG("Calling for a function evaluation");
631 time1 = tm_cpu_time();
632 #endif
633
634 /*
635 t[1]=t[0]; can't do this. lsode calls us with a different t than the t we sent in.
636 */
637 integrator_set_t(l_lsode_blsys, t[0]);
638 integrator_set_y(l_lsode_blsys, y);
639
640 #if DOTIME
641 time2 = tm_cpu_time();
642 #endif
643
644 switch(lsodesys.lastcall) {
645 case lsode_none: /* first call */
646 CONSOLE_DEBUG("FIRST CALL...");
647
648 case lsode_derivative:
649 if (lsodesys.partitioned) {
650 /* CONSOLE_DEBUG("PRE-SOLVE"); */
651 slv_presolve(l_lsode_blsys->system);
652 } else {
653 /** @TODO this doesn't ever seem to be called */
654 CONSOLE_DEBUG("RE-SOLVE");
655 slv_resolve(l_lsode_blsys->system);
656 }
657 break;
658 default:
659 case lsode_function:
660 slv_resolve(l_lsode_blsys->system);
661 break;
662 }
663
664 slv_solve(l_lsode_blsys->system);
665 slv_get_status(l_lsode_blsys->system, &status);
666 /* pass the solver status to the integrator */
667 ok = integrator_checkstatus(status);
668
669 #if DOTIME
670 time2 = tm_cpu_time() - time2;
671 #endif
672
673 if (!ok) {
674 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for derivatives");
675 /*
676 ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
677 FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");
678 for (i = 0; i< *n_eq; i++) {
679 FPRINTF(ASCERR,"y[%4d] = %f\n",i, y[i]);
680 }
681 error_reporter_end_flush();
682 */
683 lsodesys.status = lsode_nok;
684 } else {
685 lsodesys.status = lsode_ok;
686 }
687 integrator_get_ydot(l_lsode_blsys, ydot);
688
689 lsodesys.lastcall = lsode_function;
690 #if DOTIME
691 time1 = tm_cpu_time() - time1;
692 CONSOLE_DEBUG("Function evalulation has been completed in time %g. True function call time = %g",time1,time2);
693 #endif
694 }
695
696 /**
697 Evaluate the jacobian
698 */
699 static void LSODE_JEX(int *neq ,double *t, double *y,
700 int *ml ,int *mu ,double *pd, int *nrpd)
701 {
702 int nok = 0;
703 int i,j;
704
705 IntegratorLsodeData enginedata=*((IntegratorLsodeData *)l_lsode_blsys->enginedata);
706
707 UNUSED_PARAMETER(t);
708 UNUSED_PARAMETER(y);
709 UNUSED_PARAMETER(ml);
710 UNUSED_PARAMETER(mu);
711
712 /* CONSOLE_DEBUG("Calling for a gradient evaluation"); */
713 #if DOTIME
714 double time1;
715
716 CONSOLE_DEBUG("Calling for a gradient evaluation");
717 time1 = tm_cpu_time();
718 #endif
719 /*
720 * Make the real call.
721 */
722 nok = lsode_derivatives(l_lsode_blsys->system
723 , enginedata.dydx_dx
724 , enginedata.input_indices
725 , *neq
726 , enginedata.output_indices
727 , *nrpd
728 );
729
730 if (nok) {
731 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in computing the derivatives for the system. Failing...");
732 lsodesys.status = lsode_nok;
733 lsodesys.lastcall = lsode_derivative;
734 return;
735 } else {
736 lsodesys.status = lsode_ok;
737 lsodesys.lastcall = lsode_derivative;
738 }
739 /*
740 Map data from C based matrix to Fortan matrix.
741 We will send in a column major ordering vector for pd.
742 */
743 for (j=0;j<*neq;j++) { /* loop through columnns */
744 for (i=0;i<*nrpd;i++){ /* loop through rows */
745 /* CONSOLE_DEBUG("JAC[r=%d,c=%d]=%f",i,j,enginedata.dydx_dx[i][j]); */
746 *pd++ = enginedata.dydx_dx[i][j];
747 }
748 }
749
750 #if DOTIME
751 time1 = tm_cpu_time() - time1;
752 CONSOLE_DEBUG("Time to do gradient evaluation %g",time1);
753 #endif
754
755 return;
756 }
757
758 /**
759 The public function: here we do the actual integration, I guess.
760 */
761 int integrator_lsode_solve(IntegratorSystem *blsys
762 , unsigned long start_index, unsigned long finish_index
763 ){
764 slv_status_t status;
765 slv_parameters_t params;
766 IntegratorLsodeData *d;
767
768 double x[2];
769 double xend,xprev;
770 unsigned long nsamples, neq;
771 long nobs;
772 int itol, itask, mf, lrw, liw;
773 unsigned long index;
774 int istate, iopt;
775 double * rwork;
776 int * iwork;
777 double *y, *abtol, *reltol, *obs, *dydx;
778 int my_neq;
779 FILE *y_out =NULL;
780 FILE *obs_out =NULL;
781 int reporterstatus;
782
783 /* store the local variable so that we can get at stuff from inside LSODE_FEX. */
784 l_lsode_blsys = blsys;
785
786 d = (IntegratorLsodeData *)(blsys->enginedata);
787
788 /* the numer of equations must be equal to blsys->n_y, the number of states */
789 d->n_eqns = blsys->n_y;
790 assert(d->n_eqns>0);
791
792 d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
793 d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
794 d->dydx_dx = lsode_densematrix_create(d->n_eqns,d->n_eqns);
795
796 d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);
797 d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);
798
799 integrator_lsode_setup_diffs(blsys);
800
801 /* this is a lie, but we will keep it.
802 We handle any linsol/linsolqr based solver. */
803 if (strcmp(slv_solver_name(slv_get_selected_solver(blsys->system)),"QRSlv") != 0) {
804 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv must be selected before integration.");
805 return 0;
806 }
807
808 slv_get_status(l_lsode_blsys->system, &status);
809
810 if (status.struct_singular) {
811 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system is structurally singular.");
812 lsodesys.status = lsode_nok;
813 return 0;
814 }
815
816 #if defined(STATIC_LSOD) || defined (DYNAMIC_LSOD)
817
818 /* here we assume integrators.c is in charge of dynamic loading */
819
820 slv_get_parameters(blsys->system,&params);
821 lsodesys.partitioned = 1;
822
823 nsamples = integrator_getnsamples(blsys);
824 if (nsamples <2) {
825 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined.");
826 lsodesys.status = lsode_nok;
827 return 0;
828 }
829 neq = blsys->n_y;
830 nobs = blsys->n_obs;
831
832 x[0] = integrator_get_t(blsys);
833 x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */
834 lrw = 22 + 9*neq + neq*neq;
835 rwork = ASC_NEW_ARRAY_CLEAR(double, lrw+1);
836 liw = 20 + neq;
837 iwork = ASC_NEW_ARRAY_CLEAR(int, liw+1);
838 y = integrator_get_y(blsys, NULL);
839 reltol = lsode_get_rtol(blsys);
840 abtol = lsode_get_atol(blsys);
841 obs = integrator_get_observations(blsys, NULL);
842 dydx = ASC_NEW_ARRAY_CLEAR(double, neq+1);
843 if (!y || !obs || !abtol || !reltol || !rwork || !iwork || !dydx) {
844 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
845 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for lsode.");
846 lsodesys.status = lsode_nok;
847 return 0;
848 }
849
850 /*
851 Prepare args and call lsode.
852 */
853 itol = 4;
854 itask = 1;
855 istate = 1;
856 iopt = 1;
857 rwork[4] = integrator_get_stepzero(blsys);
858 rwork[5] = integrator_get_maxstep(blsys);
859 rwork[6] = integrator_get_minstep(blsys);
860 iwork[5] = integrator_get_maxsubsteps(blsys);
861 mf = 21; /* 21 = BDF with exact jacobian. 22 = BDF with finite diff Jacobian */
862
863 /* put the values from derivative system into the record */
864 integrator_setsample(blsys, start_index, x[0]);
865
866 integrator_output_init(blsys);
867
868 my_neq = (int)neq;
869
870 /*
871 First time entering lsode, x is input. After that,
872 lsode uses x as output (y output is y(x)). To drive
873 the loop ahead in time, all we need to do is keep upping
874 xend.
875 */
876 blsys->currentstep = 0;
877 for (index = start_index; index < finish_index; index++, blsys->currentstep++) {
878 xend = integrator_getsample(blsys, index+1);
879 xprev = x[0];
880 /* CONSOLE_DEBUG("BEFORE %lu LSODE CALL\n", index); */
881
882 # ifndef NO_SIGNAL_TRAPS
883 if (setjmp(g_fpe_env)==0) {
884 # endif /* NO_SIGNAL_TRAPS */
885
886 /* CONSOLE_DEBUG("Calling LSODE with end-time = %f",xend); */
887 /*
888 switch(mf){
889 case 10:
890 CONSOLE_DEBUG("Non-stiff (Adams) method; no Jacobian will be used"); break;
891 case 21:
892 CONSOLE_DEBUG("Stiff (BDF) method, user-supplied full Jacobian"); break;
893 case 22:
894 CONSOLE_DEBUG("Stiff (BDF) method, internally generated full Jacobian"); break;
895 case 24:
896 CONSOLE_DEBUG("Stiff (BDF) method, user-supplied banded jacobian"); break;
897 case 25:
898 CONSOLE_DEBUG("Stiff (BDF) method, internally generated banded jacobian"); break;
899 default:
900 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid method id %d for LSODE",mf);
901 return 0; * failure *
902 }
903 */
904
905 LSODE(&(LSODE_FEX), &my_neq, y, x, &xend,
906 &itol, reltol, abtol, &itask, &istate,
907 &iopt ,rwork, &lrw, iwork, &liw, &(LSODE_JEX), &mf);
908
909
910 # ifndef NO_SIGNAL_TRAPS
911 } else {
912 FPRINTF(stderr,
913 "Integration terminated due to float error in LSODE call.\n");
914 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
915 lsodesys.status = lsode_ok; /* clean up before we go */
916 lsodesys.lastcall = lsode_none;
917 if (y_out!=NULL) {
918 fclose(y_out);
919 }
920 if (obs_out!=NULL) {
921 fclose(obs_out);
922 }
923 return 0;
924 }
925 # endif /* NO_SIGNAL_TRAPS */
926
927 /* CONSOLE_DEBUG("AFTER %lu LSODE CALL\n", index); */
928 /* this check is better done in fex,jex, but lsode takes no status */
929 if (Solv_C_CheckHalt()) {
930 if (istate >= 0) {
931 istate=-7;
932 }
933 }
934
935 if (istate < 0 ) {
936 /* some kind of error occurred... */
937 ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
938 lsode_write_istate(istate);
939 FPRINTF(ASCERR, "\nFurthest point reached was t = %g.\n",x[0]);
940 error_reporter_end_flush();
941
942 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
943 integrator_output_close(blsys);
944 return 0;
945 }
946
947 if (lsodesys.status==lsode_nok) {
948 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to an error in derivative computations.");
949 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
950 lsodesys.status = lsode_ok; /* clean up before we go */
951 lsodesys.lastcall = lsode_none;
952 integrator_output_close(blsys);
953 return 0;
954 }
955
956 integrator_setsample(blsys, index+1, x[0]);
957 /* record when lsode actually came back */
958 integrator_set_t(blsys, x[0]);
959 integrator_set_y(blsys, y);
960 /* put x,y in d in case lsode got x,y by interpolation, as it does */
961
962 reporterstatus = integrator_output_write(blsys);
963
964 if(reporterstatus==0){
965 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration cancelled");
966 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
967 lsodesys.status = lsode_ok;
968 lsodesys.lastcall = lsode_none;
969 integrator_output_close(blsys);
970 return 0;
971 }
972
973 if (nobs > 0) {
974 # ifndef NO_SIGNAL_TRAPS
975 if (setjmp(g_fpe_env)==0) {
976 # endif /* NO_SIGNAL_TRAPS */
977
978 /* solve for obs since d isn't necessarily already
979 computed there though lsode's x and y may be.
980 Note that since lsode usually steps beyond xend
981 x1 usually wouldn't be x0 precisely if the x1/x0
982 scheme worked, which it doesn't anyway. */
983
984 LSODE_FEX(&my_neq, x, y, dydx);
985
986 /* calculate observations, if any, at returned x and y. */
987 obs = integrator_get_observations(blsys, obs);
988
989 integrator_output_write_obs(blsys);
990
991 # ifndef NO_SIGNAL_TRAPS
992 } else {
993 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in LSODE FEX call.");
994 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
995 lsodesys.status = lsode_ok; /* clean up before we go */
996 lsodesys.lastcall = lsode_none;
997 integrator_output_close(blsys);
998 return 0;
999 }
1000 # endif /* NO_SIGNAL_TRAPS */
1001 }
1002 /* CONSOLE_DEBUG("Integration completed from %3g to %3g.",xprev,x[0]); */
1003 }
1004
1005 CONSOLE_DEBUG("...");
1006 CONSOLE_DEBUG("Number of steps taken: %1d.", iwork[10]);
1007 CONSOLE_DEBUG("Number of function evaluations: %1d.", iwork[11]);
1008 CONSOLE_DEBUG("Number of Jacobian evaluations: %1d.", iwork[12]);
1009 CONSOLE_DEBUG("...");
1010
1011
1012 lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1013
1014 /*
1015 * return the system to its original state.
1016 */
1017
1018 lsodesys.status = lsode_ok;
1019 lsodesys.lastcall = lsode_none;
1020
1021 integrator_output_close(blsys);
1022
1023 CONSOLE_DEBUG("--- LSODE done ---");
1024 return 1;
1025
1026 #else /* STATIC_LSOD || DYNAMIC_LSOD */
1027
1028 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration will not be performed. LSODE binary not available.");
1029 lsodesys.status = lsode_nok;
1030 return 0;
1031
1032 #endif
1033 }
1034
1035 /**
1036 Function XASCWV is an error reporting function replacing the XERRWV
1037 routine in lsode.f. The call signature is the same with the original Fortran
1038 function.
1039
1040 @see the comments for 'xerrwv' from lsode.f, with which XASCWV is compatible...
1041
1042 @param msg = the message (hollerith literal or integer array).
1043 @param nmes = the length of msg (number of characters).
1044 @param nerr = the error number (not used).
1045 @param level = the error level..
1046 0 or 1 means recoverable (control returns to caller).
1047 2 means fatal (run is aborted--see note below).
1048 @param ni = number of integers (0, 1, or 2) to be printed with message.
1049 @param i1,i2 = integers to be printed, depending on ni.
1050 @param nr = number of reals (0, 1, or 2) to be printed with message.
1051 @param r1,r2 = reals to be printed, depending on nr.
1052 */
1053 void XASCWV( char *msg, /* pointer to start of message */
1054 int *nmes, /* the length of msg (number of characters) */
1055 int *nerr, /* the error number (not used). */
1056 int *level,
1057 int *ni,
1058 int *i1,
1059 int *i2,
1060 int *nr,
1061 double *r1,
1062 double *r2
1063 ){
1064 switch(*nerr){
1065 case 204:
1066 if(*nr==0)return;
1067 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);
1068 break;
1069 case 205:
1070 if(*nr==0)return;
1071 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);
1072 break;
1073
1074 default:
1075 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1076
1077 /* note that %.*s means that a string length (integer) and string pointer are being required */
1078 FPRINTF(stderr,"LSODE error: %.*s",*nmes,msg);
1079 if (*ni == 1) {
1080 FPRINTF(stderr,"\nwhere i1 = %d\n",*i1);
1081 }
1082 if (*ni == 2) {
1083 FPRINTF(stderr,"\nwhere i1 = %d, i2 = %d",*i1,*i2);
1084 }
1085 if (*nr == 1) {
1086 FPRINTF(stderr,"\nwhere r1 = %.13g", *r1);
1087 }
1088 if (*nr == 2) {
1089 FPRINTF(stderr,"\nwhere r1 = %.13g, r2 = %.13g", *r1,*r2);
1090 }
1091 error_reporter_end_flush();
1092 }
1093
1094 if (*level != 2) {
1095 return;
1096 }
1097
1098 /* NOT reached. lsode does NOT make level 2 calls in our version. */
1099 error_reporter_end_flush();
1100 Asc_Panic(3,"xascwv", "LSODE really really confused");
1101 }

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