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

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