/[ascend]/trunk/solvers/lsode/asc_lsode.c
ViewVC logotype

Contents of /trunk/solvers/lsode/asc_lsode.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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