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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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