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

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