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

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