/[ascend]/trunk/tcltk98/generic/interface/Lsode.c
ViewVC logotype

Contents of /trunk/tcltk98/generic/interface/Lsode.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 428 - (show annotations) (download) (as text)
Tue Apr 4 06:41:25 2006 UTC (16 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 24816 byte(s)
Much progress on the MinGW build of the Tcl/Tk interface. Builds and links now,
just having some problems getting the DLLs for Tcl/Tk to resolve at runtime.
1 /*
2 * Lsode.c
3 * by Kirk Abbott and Ben Allan
4 * Created: 1/94
5 * Version: $Revision: 1.29 $
6 * Version control file: $RCSfile: Lsode.c,v $
7 * Date last modified: $Date: 2000/01/25 02:26:31 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the ASCEND Tcl/Tk interface
11 *
12 * Copyright 1997, Carnegie Mellon University
13 *
14 * The ASCEND Tcl/Tk interface is free software; you can redistribute
15 * it and/or modify it under the terms of the GNU General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * The ASCEND Tcl/Tk interface is distributed in hope that it will be
20 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with the program; if not, write to the Free Software Foundation,
26 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
27 * COPYING. COPYING is found in ../compiler.
28 */
29
30 #ifndef NO_SIGNAL_TRAPS
31 #include <signal.h>
32 #include <setjmp.h>
33 #endif /* NO_SIGNAL_TRAPS */
34 #include <tcl.h> /* for Integrators.h, Sensitivity.h*/
35 #include <utilities/ascConfig.h>
36 #include <utilities/ascSignal.h>
37 #include <utilities/ascMalloc.h>
38 #include <utilities/ascPanic.h>
39 #include <general/tm_time.h>
40 #include <general/list.h>
41 #include <compiler/compiler.h>
42 #include <compiler/instance_enum.h>
43 #include <compiler/symtab.h>
44 #include <compiler/fractions.h>
45 #include <compiler/dimen.h>
46 #include <compiler/parentchild.h>
47 #include <compiler/module.h>
48 #include <compiler/library.h>
49 #include <compiler/child.h>
50 #include <compiler/type_desc.h>
51 #include <compiler/instance_name.h>
52 #include <compiler/atomvalue.h>
53 #include <solver/slv_types.h>
54 #include <solver/mtx.h>
55 #include <solver/var.h>
56 #include <solver/rel.h>
57 #include <solver/discrete.h>
58 #include <solver/conditional.h>
59 #include <solver/logrel.h>
60 #include <solver/bnd.h>
61 #include <solver/slv_common.h>
62 #include <solver/linsol.h>
63 #include <solver/linsolqr.h>
64 #include <solver/slv_client.h>
65 #include <compiler/types.h>
66 #include <compiler/functype.h>
67 #include <compiler/func.h>
68 #include <compiler/extfunc.h>
69 #include <compiler/extcall.h>
70 #include <compiler/relation_type.h>
71 #include "old_utils.h"
72 #include "Integrators.h"
73 #include "Lsode.h"
74 #include "Sensitivity.h" /* see the packages dir */
75
76 #ifndef lint
77 static CONST char LsodeID[] = "$Id: Lsode.c,v 1.29 2000/01/25 02:26:31 ballan Exp $";
78 #endif
79
80
81 /*
82 * NOUNDERBARS --> FORTRAN compiler naming convention for subroutine
83 * is wierd. WIN32/CRAY is treated as special case
84 */
85 #ifdef APOLLO
86 #define NOUNDERBARS TRUE
87 #endif
88 #ifdef _HPUX_SOURCE
89 #define NOUNDERBARS TRUE
90 #endif
91 /* AIX xlf will not suffix an underbar on a symbol
92 * unless xlf is given the ``-qextname'' option
93 */
94 #ifdef _AIX
95 #define NOUNDERBARS TRUE
96 #endif
97
98 #ifdef NOUNDERBARS
99 #define LSODE lsode
100 #define LSODE_JEX jex
101 #define LSODE_FEX fex
102 #define GETCOMMON get_lsode_common
103 #define XASCWV xascwv
104 #else
105 /* sun, __alpha, __sgi, ... */
106 #define LSODE lsode_
107 #define LSODE_JEX jex_
108 #define LSODE_FEX fex_
109 #define GETCOMMON get_lsode_common_
110 #define XASCWV xascwv_
111 #endif
112
113 #if defined(CRAY) || (defined(__WIN32__) && !defined(__MINGW32_VERSION))
114 #undef LSODE
115 #undef LSODE_JEX
116 #undef LSODE_FEX
117 #undef GETCOMMON
118 #undef XASCWV
119 #define XASCWV XASCWV
120 #define LSODE LSODE
121 #define LSODE_JEX JEX
122 #define LSODE_FEX FEX
123 #define GETCOMMON GET_LSODE_COMMON
124 #endif
125
126 #define DOTIME FALSE
127
128 /* definitions of lsode supported children of atoms, etc */
129 /********************************************************************/
130 /* default input tolerances for lsode */
131 #define RTOLDEF 1e-6
132 #define ATOLDEF 1e-6
133 /* solver_var children expected for state variables */
134 static symchar *g_symbols[2];
135 #define STATERTOL g_symbols[0]
136 #define STATEATOL g_symbols[1]
137 static
138 void InitTolNames(void)
139 {
140 STATERTOL = AddSymbol("ode_rtol");
141 STATEATOL = AddSymbol("ode_atol");
142 }
143
144 /***** interface implementation notes *****/
145 /*
146 * As fortran io is unreliably portable (vc5+digital fortran)
147 * we have converted xerrwv to xascwv provided here.
148 *
149 * The lsode interface variable t is actually an array of
150 * 2 doubles rather than just 1. The first is the the one
151 * used by lsode. The second is used by LSODE_FEX to tell
152 * what the last time it was called was. This is so the
153 * C driver can tell if it needs to resolve d to compute
154 * observation variables. If x[0]==x[1] we save ourselves
155 * a solve.
156 ! Note!!!!
157 ! The above doesn't work since lsode doesn't use the same
158 ! t internally that we hand it.
159 */
160
161 enum Lsode_enum {
162 lsode_none, /* true on first call */
163 lsode_function, lsode_derivative, /* functions or gradients done */
164 lsode_sparse, lsode_dense, /* what type of backend should we */
165 lsode_band, /* use for the integrator */
166 lsode_ok, lsode_nok /* bad return from func or grad */
167 };
168
169 static struct Lsode_Data {
170 enum Lsode_enum lastcall; /* type of last call; func or grad */
171 enum Lsode_enum status; /* solve status */
172 int partitioned; /* partioned func evals or not */
173 struct rel_relation **rlist; /* relations to differentiate */
174 slv_system_t sys; /* the main solve system */
175 } lsodesys = {lsode_none, lsode_ok, 1, NULL, NULL};
176
177
178 /*
179 * Type of function used to evaluate derivative system.
180 */
181 typedef void (*fex_t)( int * ,double * ,double *,double *);
182
183 /*
184 * Type of function used to evaluate jacobian system.
185 */
186 typedef void (*jex_t)(int *, double *, double *, int *, int *,
187 double *, int *);
188
189 /*
190 void LSODE(&fex, &neq, y, &x, &xend, &itol, reltol, abtol, &itask,
191 &istate, &iopt ,rwork, &lrw, iwork, &liw, &jex, &mf);
192
193 This is the FORTRAN call to LSODE.
194
195 Removed 'extern' so that linker will complain if it's not there.
196 */
197 void LSODE( fex_t
198 ,int *
199 ,double *
200 ,double *
201 ,double *
202 ,int *
203 ,double *
204 ,double *
205 ,int *
206 ,int *
207 ,int *
208 ,double *
209 ,int *
210 ,int *
211 ,int *
212 ,jex_t
213 ,int *
214 );
215
216 /********************************************************************/
217 /* allocates, fills, and returns the atol vector based on BLSODE */
218 /* State variables missing child ode_rtol will be defaulted to ATOLDEF */
219
220 static double *blsode_get_atol( struct Integ_system_t *blsys) {
221
222 struct Instance *tol;
223 double *atoli;
224 int i,len;
225
226 len = blsys->n_y;
227 atoli = (double *)ascmalloc((blsys->n_y+1)*sizeof(double));
228 if (atoli == NULL) {
229 FPRINTF(stderr,"ERROR: (blsode_get_atol) Insufficient memory.\n");
230 return atoli;
231 }
232 InitTolNames();
233 for (i=0; i<len; i++) {
234 tol = ChildByChar(var_instance(blsys->y[i]),STATEATOL);
235 if (tol == NULL || !AtomAssigned(tol) ) {
236 atoli[i] = ATOLDEF;
237 FPRINTF(stderr,"WARNING: (blsode_get_atol) Assuming atol = %3g\n",
238 ATOLDEF);
239 FPRINTF(stderr, "for ode_atol child undefined for state variable %ld.\n",
240 blsys->y_id[i]);
241 } else {
242 atoli[i] = RealAtomValue(tol);
243 }
244 }
245 atoli[len] = ATOLDEF;
246 return atoli;
247 }
248
249 /********************************************************************/
250 /* Allocates, fills, and returns the rtol vector based on BLSODE */
251 /* State variables missing child ode_rtol will be defaulted to RTOLDEF */
252 static double *blsode_get_rtol( struct Integ_system_t *blsys) {
253
254 struct Instance *tol;
255 double *rtoli;
256 int i,len;
257
258 len = blsys->n_y;
259 rtoli = (double *)ascmalloc((blsys->n_y+1)*sizeof(double));
260 if (rtoli == NULL) {
261 FPRINTF(stderr,"ERROR: (blsode_get_rtol) Insufficient memory.\n");
262 return rtoli;
263 }
264 InitTolNames();
265 for (i=0; i<len; i++) {
266 tol = ChildByChar(var_instance(blsys->y[i]),STATERTOL);
267 if (tol == NULL || !AtomAssigned(tol) ) {
268 rtoli[i] = RTOLDEF;
269 FPRINTF(stderr,"WARNING: (blsode_get_rtol) Assuming rtol = %3g\n",
270 ATOLDEF);
271 FPRINTF(stderr, "for ode_rtol child undefined for state variable %ld.\n",
272 blsys->y_id[i]);
273 } else {
274 rtoli[i] = RealAtomValue(tol);
275 }
276 }
277 rtoli[len] = RTOLDEF;
278 return rtoli;
279 }
280
281 /********************************************************************/
282
283 static void write_istate( int istate) {
284 FPRINTF(stderr,"lsode integrator:");
285 switch (istate) {
286 case -1:
287 FPRINTF(stderr,"Excess steps taken on this call (perhaps wrong MF).\n");
288 break;
289 case -2:
290 FPRINTF(stderr,"Excess accuracy requested (tolerances too small).\n");
291 break;
292 case -3:
293 FPRINTF(stderr,"Illegal input detected.\n");
294 break;
295 case -4:
296 FPRINTF(stderr,"Repeated error test failures.\n");
297 break;
298 case -5:
299 FPRINTF(stderr,"Repeated convergence failures.\n");
300 break;
301 case -6:
302 FPRINTF(stderr,"Error weight became zero during problem.\n");
303 break;
304 case -7:
305 FPRINTF(stderr,"User patience became zero during problem.\n");
306 break;
307 default:
308 FPRINTF(stderr,"Unknown error code %d\n",istate);
309 break;
310 }
311 }
312 /********************************************************************/
313
314 static void freeMem(double *y, double *reltol, double *abtol, double *rwork,
315 int *iwork, double *obs, double *dydx)
316 {
317 if (y != NULL) {
318 ascfree((double *)y);
319 }
320 if (reltol != NULL) {
321 ascfree((double *)reltol);
322 }
323 if (abtol != NULL) {
324 ascfree((double *)abtol);
325 }
326 if (rwork != NULL) {
327 ascfree((double *)rwork);
328 }
329 if (iwork != NULL) {
330 ascfree((int *)iwork);
331 }
332 if (obs != NULL) {
333 ascfree((double *)obs);
334 }
335 if (dydx != NULL) {
336 ascfree((double *)dydx);
337 }
338 }
339
340 /********************************************************************/
341 /*
342 * The current way that we are getting the derivatives (if the problem
343 * was solved partitioned) messes up the slv_system so that we *have*
344 * to do a *presolve* rather than a simply a *resolve* before doing
345 * function calls. This code below attempts to handle these cases.
346 */
347 static void LSODE_FEX( int *n_eq ,double *t ,double *y ,double *ydot)
348 {
349 slv_status_t status;
350 /* slv_parameters_t parameters; pity lsode doesn't allow error returns */
351 int i;
352 unsigned long ok;
353
354 #if DOTIME
355 double time1,time2;
356 #endif
357
358 #if DOTIME
359 FPRINTF(stderr,"\n** Calling for a function evaluation **\n");
360 time1 = tm_cpu_time();
361 #endif
362
363 /* t[1]=t[0]; can't do this. lsode calls us with a different t
364 * than the x we sent in.
365 */
366 Asc_IntegSetDX(t[0]);
367 Asc_IntegSetDY(y);
368
369 #if DOTIME
370 time2 = tm_cpu_time();
371 #endif
372
373 switch(lsodesys.lastcall) {
374 case lsode_none: /* first call */
375 case lsode_derivative:
376 if (lsodesys.partitioned) {
377 slv_presolve(lsodesys.sys);
378 } else {
379 slv_resolve(lsodesys.sys);
380 }
381 break;
382 default:
383 case lsode_function:
384 slv_resolve(lsodesys.sys);
385 break;
386 }
387 slv_solve(lsodesys.sys);
388 slv_get_status(lsodesys.sys, &status);
389 ok = Asc_IntegCheckStatus(status);
390
391 #if DOTIME
392 time2 = tm_cpu_time() - time2;
393 #endif
394
395 if (!ok) {
396 FPRINTF(stderr,"Unable to compute the vector of derivatives with the\n");
397 FPRINTF(stderr,"following values for the state variables:\n");
398 for (i = 0; i< *n_eq; i++) {
399 FPRINTF(stderr,"y[%4d] = %f\n",i, y[i]);
400 }
401 FPRINTF(stderr,"\n");
402 lsodesys.status = lsode_nok;
403 } else {
404 lsodesys.status = lsode_ok;
405 }
406 Asc_IntegGetDDydx(ydot);
407
408 lsodesys.lastcall = lsode_function;
409 #if DOTIME
410 time1 = tm_cpu_time() - time1;
411 FPRINTF(stderr,"** Function evalulation has been completed in %g**\n",time1);
412 FPRINTF(stderr,"** True function call time = %g**\n",time2);
413 #endif
414 }
415
416 /********************************************************************/
417
418 static void LSODE_JEX(int *neq ,double *t, double *y,
419 int *ml ,int *mu ,double *pd, int *nrpd)
420 {
421 int nok = 0;
422 int i,j;
423
424 (void)t; /* stop gcc whine about unused parameter */
425 (void)y; /* stop gcc whine about unused parameter */
426 (void)ml; /* stop gcc whine about unused parameter */
427 (void)mu; /* stop gcc whine about unused parameter */
428
429 #if DOTIME
430 double time1;
431 #endif
432
433 #if DOTIME
434 FPRINTF(stderr,"\n** Calling for a gradient evaluation **\n");
435 time1 = tm_cpu_time();
436 #endif
437 /*
438 * Make the real call.
439 */
440 nok = Asc_BLsodeDerivatives(lsodesys.sys, g_intg_diff.dydx_dx,
441 g_intg_diff.input_indices, *neq,
442 g_intg_diff.output_indices, *nrpd);
443 if (nok) {
444 FPRINTF(stderr,"Error in computing the derivatives for the system\n");
445 FPRINTF(stderr,"Failing...\n");
446 lsodesys.status = lsode_nok;
447 lsodesys.lastcall = lsode_derivative;
448 return;
449 } else {
450 lsodesys.status = lsode_ok;
451 lsodesys.lastcall = lsode_derivative;
452 }
453 /*
454 * Map data from C based matrix to Fortan matrix.
455 * We will send in a column major ordering vector for pd.
456 */
457 for (j=0;j<*neq;j++) {
458 for (i=0;i<*nrpd;i++) { /* column wise - hence rows change faster */
459 *pd++ = g_intg_diff.dydx_dx[i][j];
460 }
461 }
462
463 #if DOTIME
464 time1 = tm_cpu_time() - time1;
465 FPRINTF(stderr,"********** Time to do gradient evaluation %g\n",time1);
466 #endif
467
468 return;
469 }
470
471 /********************************************************************/
472
473
474 void Asc_BLsodeIntegrate(slv_system_t sys, unsigned long start_index ,
475 unsigned long finish_index, struct Integ_system_t *blsys)
476 {
477 slv_status_t status;
478 slv_parameters_t params;
479 double x[2];
480 double xend,xprev;
481 unsigned long nsamples, neq;
482 long nobs;
483 int itol, itask, mf, lrw, liw;
484 unsigned long index;
485 int istate, iopt;
486 double * rwork;
487 int * iwork;
488 double *y, *abtol, *reltol, *obs, *dydx;
489 int my_neq;
490 FILE *y_out =NULL;
491 FILE *obs_out =NULL;
492
493 lsodesys.sys = sys;
494 slv_get_status(lsodesys.sys, &status);
495
496 if (status.struct_singular) {
497 FPRINTF(stderr, "\n");
498 FPRINTF(stderr, "Integration will not be performed.\n");
499 FPRINTF(stderr, "The system is structurally singular.\n");
500 lsodesys.status = lsode_nok;
501 return;
502 }
503
504 #if defined(STATIC_LSOD) || defined (DYNAMIC_LSOD)
505 /* here we assume integrators.c is in charge of dynamic loading */
506
507 slv_get_parameters(lsodesys.sys,&params);
508 lsodesys.partitioned = 1;
509
510 nsamples = Asc_IntegGetNSamples();
511 if (nsamples <2) {
512 FPRINTF(stderr, "\n");
513 FPRINTF(stderr, "Integration will not be performed.\n");
514 FPRINTF(stderr, "The system has no end sample time defined.\n");
515 lsodesys.status = lsode_nok;
516 return;
517 }
518 neq = blsys->n_y;
519 nobs = blsys->n_obs;
520
521 x[0] = Asc_IntegGetDX();
522 x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */
523 lrw = 22 + 9*neq + neq*neq;
524 rwork = (double *)asccalloc(lrw+1, sizeof(double));
525 liw = 20 + neq;
526 iwork = (int *)asccalloc(liw+1, sizeof(int));
527 y = Asc_IntegGetDY(NULL);
528 reltol = blsode_get_rtol(blsys);
529 abtol = blsode_get_atol(blsys);
530 obs = Asc_IntegGetDObs(NULL);
531 dydx = (double *)asccalloc(neq+1, sizeof(double));
532 if (!y || !obs || !abtol || !reltol || !rwork || !iwork || !dydx) {
533 freeMem(y,reltol,abtol,rwork,iwork,obs,dydx);
534 FPRINTF(stderr,"Insufficient memory for blsode.\n");
535 lsodesys.status = lsode_nok;
536 return;
537 }
538
539 y_out = Asc_IntegOpenYFile();
540 obs_out = Asc_IntegOpenObsFile();
541 Asc_IntegReleaseYFile();
542 Asc_IntegReleaseObsFile();
543 /*
544 * Prepare args and call lsode.
545 */
546 itol = 4;
547 itask = 1;
548 istate = 1;
549 iopt = 1;
550 rwork[4] = Asc_IntegGetStepZero(blsys);
551 rwork[5] = Asc_IntegGetStepMax(blsys);
552 rwork[6] = Asc_IntegGetStepMin(blsys);
553 iwork[5] = Asc_IntegGetMaxSteps(blsys);
554 mf = 21; /* gradients 22 -- implies finite diff */
555
556 /* put the values from derivative system into the record */
557 Asc_IntegSetXSamplei(start_index, x[0]);
558 /* write headers to yout, obsout and initial points */
559 Asc_IntegPrintYHeader(y_out,blsys);
560 Asc_IntegPrintYLine(y_out,blsys);
561 Asc_IntegPrintObsHeader(obs_out,blsys);
562 Asc_IntegPrintObsLine(obs_out,blsys);
563
564 my_neq = (int)neq;
565 /*
566 * First time entering lsode, x is input. After that,
567 * lsode uses x as output (y output is y(x)). To drive
568 * the loop ahead in time, all we need to do is keep upping
569 * xend.
570 */
571 for (index = start_index; index < finish_index; index++) {
572 xend = Asc_IntegGetXSamplei(index+1);
573 xprev = x[0];
574 print_debug("BEFORE %lu BLSODE CALL\n", index);
575
576 #ifndef NO_SIGNAL_TRAPS
577 if (setjmp(g_fpe_env)==0) {
578 #endif /* NO_SIGNAL_TRAPS */
579 LSODE(&(LSODE_FEX), &my_neq, y, x, &xend,
580 &itol, reltol, abtol, &itask, &istate,
581 &iopt ,rwork, &lrw, iwork, &liw, &(LSODE_JEX), &mf);
582 #ifndef NO_SIGNAL_TRAPS
583 } else {
584 FPRINTF(stderr,
585 "Integration terminated due to float error in LSODE call.\n");
586 freeMem(y,reltol,abtol,rwork,iwork,obs,dydx);
587 lsodesys.status = lsode_ok; /* clean up before we go */
588 lsodesys.lastcall = lsode_none;
589 if (y_out!=NULL) {
590 fclose(y_out);
591 }
592 if (obs_out!=NULL) {
593 fclose(obs_out);
594 }
595 return;
596 }
597 #endif /* NO_SIGNAL_TRAPS */
598
599 print_debug("AFTER %lu LSODE CALL\n", index);
600 /* this check is better done in fex,jex, but lsode takes no status */
601 if (Solv_C_CheckHalt()) {
602 if (istate >= 0) {
603 istate=-7;
604 }
605 }
606 if (istate < 0 ) {
607 FPRINTF(stderr, "!! Asc_BLsodeIntegrate: ");
608 FPRINTF(stderr, "index = %lu ", index);
609 FPRINTF(stderr, "istate = %d ", istate);
610 FPRINTF(stderr, "farthest point reached: x = %g\n",x[0]);
611 write_istate(istate);
612 freeMem(y,reltol,abtol,rwork,iwork,obs,dydx);
613 if (y_out!=NULL) {
614 fclose(y_out);
615 }
616 if (obs_out!=NULL) {
617 fclose(obs_out);
618 }
619 return;
620 }
621
622 if (lsodesys.status==lsode_nok) {
623 FPRINTF(stderr,
624 "Integration terminated due to an error in derivative computations.\n"
625 );
626 freeMem(y,reltol,abtol,rwork,iwork,obs,dydx);
627 lsodesys.status = lsode_ok; /* clean up before we go */
628 lsodesys.lastcall = lsode_none;
629 if (y_out!=NULL) {
630 fclose(y_out);
631 }
632 if (obs_out!=NULL) {
633 fclose(obs_out);
634 }
635 return;
636 }
637
638 Asc_IntegSetXSamplei(index+1, x[0]);
639 /* record when lsode actually came back */
640 Asc_IntegSetDX(x[0]);
641 Asc_IntegSetDY(y);
642 /* put x,y in d in case lsode got x,y by interpolation, as it does */
643 Asc_IntegPrintYLine(y_out,blsys);
644 if (nobs > 0) {
645 #ifndef NO_SIGNAL_TRAPS
646 if (setjmp(g_fpe_env)==0) {
647 #endif /* NO_SIGNAL_TRAPS */
648 /* solve for obs since d isn't necessarily already
649 computed there though lsode's x and y may be.
650 Note that since lsode usually steps beyond xend
651 x1 usually wouldn't be x0 precisely if the x1/x0
652 scheme worked, which it doesn't anyway. */
653 LSODE_FEX(&my_neq, x, y, dydx);
654 /* calculate observations, if any, at returned x and y. */
655 obs = Asc_IntegGetDObs(obs);
656 Asc_IntegPrintObsLine(obs_out,blsys);
657 #ifndef NO_SIGNAL_TRAPS
658 } else {
659 FPRINTF(stderr,
660 "Integration terminated due to float error in LSODE FEX call.\n"
661 );
662 freeMem(y,reltol,abtol,rwork,iwork,obs,dydx);
663 lsodesys.status = lsode_ok; /* clean up before we go */
664 lsodesys.lastcall = lsode_none;
665 if (y_out!=NULL) {
666 fclose(y_out);
667 }
668 if (obs_out!=NULL) {
669 fclose(obs_out);
670 }
671 return;
672 }
673 #endif /* NO_SIGNAL_TRAPS */
674 }
675 FPRINTF(stdout, "Integration completed from ");
676 FPRINTF(stdout, "%3g to %3g\n",xprev,x[0]);
677 }
678
679 FPRINTF(stdout, "\n");
680 FPRINTF(stdout, "Number of steps taken: %1d\n", iwork[10]);
681 FPRINTF(stdout, "Number of function evaluations: %1d\n", iwork[11]);
682 FPRINTF(stdout, "Number of Jacobian evaluations: %1d\n", iwork[12]);
683 FPRINTF(stdout, "\n");
684
685 freeMem(y,reltol,abtol,rwork,iwork,obs,dydx);
686 /*
687 * return the system to its original state.
688 */
689 lsodesys.status = lsode_ok;
690 lsodesys.lastcall = lsode_none;
691 if (y_out!=NULL) {
692 fclose(y_out);
693 }
694 if (obs_out!=NULL) {
695 fclose(obs_out);
696 }
697 FPRINTF(stdout, "blsode done.\n");
698
699 #else /* STATIC_LSOD || DYNAMIC_LSOD */
700
701 FPRINTF(stderr, "\n");
702 FPRINTF(stderr, "Integration will not be performed.\n");
703 FPRINTF(stderr, "LSODE binary not available.\n");
704 lsodesys.status = lsode_nok;
705 return;
706
707 #endif
708 }
709
710 /*
711 * This function is xascwv, a compatible
712 * replacement for xerrwv in lsode.f by Hindemarsh.
713 */
714 /* Ported to C, approximately, by Benjamin Allan, 9/20/97.
715 * This C source placed in the public domain, since lsode
716 * itself is PD. This placedment in no way affects the GNU
717 * licensing OF the rest OF the ascend system source code.
718 * FORTRAN header:
719 *-----------------------------------------------------------------------
720 * subroutine xerrwv, constitute
721 * a simplified version of the slatec error handling package.
722 * written by a. c. hindmarsh at llnl. version of march 30, 1987.
723 * this version is in double precision.
724 *
725 * all arguments are input arguments.
726 *
727 * msg = the message (hollerith literal or integer array).
728 * nmes = the length of msg (number of characters).
729 * nerr = the error number (not used).
730 * level = the error level..
731 * 0 or 1 means recoverable (control returns to caller).
732 * 2 means fatal (run is aborted--see note below).
733 * ni = number of integers (0, 1, or 2) to be printed with message.
734 * i1,i2 = integers to be printed, depending on ni.
735 * nr = number of reals (0, 1, or 2) to be printed with message.
736 * r1,r2 = reals to be printed, depending on nr.
737 *
738 * note.. this routine is machine-dependent and specialized for use
739 * in limited context, in the following ways..
740 * 1. the number of hollerith characters stored per word, denoted
741 * by ncpw below, is a data-loaded constant.
742 * 2. the value of nmes is assumed to be at most 60.
743 * (multi-line messages are generated by repeated calls.)
744 * 3. if level = 2, control passes to the statement stop
745 * to abort the run. this statement may be machine-dependent.
746 * 4. r1 and r2 are assumed to be in double precision and are printed
747 * in d21.13 format.
748 * 5. the common block /eh0001/ below is data-loaded (a machine-
749 * dependent feature) with default values.
750 * this block is needed for proper retention of parameters used by
751 * this routine which the user can reset by calling xsetf or xsetun.
752 * the variables in this block are as follows..
753 * mesflg = print control flag..
754 * 1 means print all messages (the default).
755 * 0 means no printing.
756 * lunit = logical unit number for messages.
757 * the default is 6 (machine-dependent).
758 *-----------------------------------------------------------------------
759 * the following are instructions for installing this routine
760 * in different machine environments.
761 *
762 * to change the default output unit, change the data statement
763 * in the block data subprogram below.
764 *
765 * for a different number of characters per word, change the
766 * data statement setting ncpw below, and format 10. alternatives for
767 * various computers are shown in comment cards.
768 *
769 * for a different run-abort command, change the statement following
770 * statement 100 at the end.
771 *-----------------------------------------------------------------------
772 */
773 void XASCWV( char *msg, /* array of char/int not NULL ended, len *nmes */
774 int *nmes, /* len of msg */
775 int *nerr,
776 int *level,
777 int *ni,
778 int *i1,
779 int *i2,
780 int *nr,
781 double *r1,
782 double *r2
783 )
784 {
785 (void)nerr;
786 /* ignore
787 * integer i,lun, lunit, mesflg, ncpw, nch, nwds
788 * common /eh0001/ mesflg, lunit
789 * data ncpw/4/
790 *
791 * ncpw = sizeof(int);
792 * if (mesflg .eq. 0) go to 100 !ignore io suppresion
793 * lun = lunit
794 *
795 * nch = min0(nmes,60)
796 * nwds = nch/ncpw
797 * if (nch .ne. nwds*ncpw) nwds = nwds + 1
798 */
799 /* write the message. lot easier in C, geez */
800 FPRINTF(stderr,"%.*s\n",*nmes,msg);
801 if (*ni == 1) {
802 FPRINTF(stderr," in above message, i1 = %d\n",*i1);
803 }
804 if (*ni == 2) {
805 FPRINTF(stderr," in above message, i1 = %d i2 = %d\n",*i1,*i2);
806 }
807 if (*nr == 1) {
808 FPRINTF(stderr," in above message, r1 = %21.13g\n", *r1);
809 }
810 if (*nr == 2) {
811 FPRINTF(stderr," above, r1 = %21.13g r2 = %21.13g\n", *r1,*r2);
812 }
813 if (*level != 2) {
814 return;
815 }
816 /* NOT reached. lsode does NOT make level 2 calls in our version. */
817 Asc_Panic(3,"xascwv", "LSODE really really confused");
818 return; /* not reached */
819 }

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