/[ascend]/trunk/tcltk/generic/interface/Integrators.c
ViewVC logotype

Contents of /trunk/tcltk/generic/interface/Integrators.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 571 - (show annotations) (download) (as text)
Tue May 9 00:14:59 2006 UTC (14 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 53278 byte(s)
Renaming 'tcltk98' to 'tcltk', continued...
1 /*
2 * Integrators.c
3 * by Kirk Abbott and Ben Allan
4 * Created: 1/94
5 * Version: $Revision: 1.32 $
6 * Version control file: $RCSfile: Integrators.c,v $
7 * Date last modified: $Date: 2003/08/23 18:43:06 $
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 /*
31 * This module defines the general integration auxillaries for Ascend.
32 */
33
34 #include <time.h>
35 #include <tcl.h>
36 #include <utilities/ascConfig.h>
37 #include <utilities/ascMalloc.h>
38 #include <utilities/readln.h>
39 #include <general/list.h>
40 #include <general/dstring.h>
41 #include <compiler/compiler.h>
42 #include <compiler/instance_enum.h>
43 #include <compiler/fractions.h>
44 #include <compiler/dimen.h>
45 #include <compiler/units.h>
46 #include <compiler/module.h>
47 #include <compiler/library.h>
48 #include <compiler/types.h>
49 #include <compiler/child.h>
50 #include <compiler/type_desc.h>
51 #include <compiler/atomvalue.h>
52 #include <compiler/instance_name.h>
53 #include <compiler/instquery.h>
54 #include <compiler/parentchild.h>
55 #include <compiler/symtab.h>
56 #include <compiler/instance_io.h>
57 #include <solver/slv_types.h>
58 #include <solver/mtx.h>
59 #include <solver/var.h>
60 /*
61 * note: the analytic jacobian routines (state matrix) depend on the
62 * assumption that struct var_variable *<--> struct Instance *.
63 */
64 #include <solver/rel.h>
65 #include <solver/discrete.h>
66 #include <solver/conditional.h>
67 #include <solver/logrel.h>
68 #include <solver/bnd.h>
69 #include <solver/slv_common.h>
70 #include <solver/linsol.h>
71 #include <solver/linsolqr.h>
72 #include <solver/slv_client.h>
73 #include "HelpProc.h"
74 #include "Integrators.h"
75 #include "BrowserQuery.h"
76 #include "Qlfdid.h"
77 #include "UnitsProc.h"
78 #include "BrowserProc.h"
79 #include "HelpProc.h"
80 #include "SolverGlobals.h"
81 #include "Lsode.h"
82
83 #ifndef lint
84 static CONST char IntegratorsID[] = "$Id: Integrators.c,v 1.32 2003/08/23 18:43:06 ballan Exp $";
85 #endif
86
87
88 #define SNULL (char *)NULL
89
90 static symchar *g_symbols[3];
91
92 /* The following names are of solver_var children or attributes
93 * we support (at least temporarily) to determine who is a state and
94 * who matching derivative.
95 * These should be supported directly in a future solveratominst.
96 */
97 #define STATEFLAG g_symbols[0]
98 /* Integer child. 0= algebraic, 1 = state, 2 = derivative, 3 = 2nd deriv etc */
99 #define STATEINDEX g_symbols[1]
100 /* Integer child. all variables with the same STATEINDEX value are taken to
101 * be derivatives of the same state variable. We really need a compiler
102 * that maintains this info by backpointers, but oh well until then.
103 */
104 #define OBSINDEX g_symbols[2]
105 /* Integer child. All variables with OBSINDEX !=0 will be recorded in
106 * the blsode output file. Tis someone else's job to grok this output.
107 */
108
109 #define INTDEBUG 0
110
111 struct Integ_var_t {
112 long index;
113 long type;
114 struct var_variable *i;
115 };
116 /* temp catcher of dynamic variable and observation variable data */
117
118 struct Integ_samples_t {
119 long ns;
120 double *x;
121 dim_type *d;
122 };
123 /* an array of sample 'times' for reference during integration */
124
125 /* global to the world */
126 struct Integ_DiffData g_intg_diff = {0L,NULL,NULL,NULL,NULL,NULL,NULL};
127
128 /*Define items for integration interface. they are global to this file.
129 this globalness is bad. */
130 static slv_system_t g_intgsys_cur = NULL;
131 /* derivative system children */
132 static struct var_variable *g_intginst_d_x = NULL;
133 /* derivative system sizes */
134 static long g_intginst_d_n_eq, g_intginst_d_n_obs;
135 /* derivative system array var lists. Not null terminated */
136 static struct var_variable **g_intginst_d_y_vars = NULL;
137 static struct var_variable **g_intginst_d_dydx_vars = NULL;
138 static struct var_variable **g_intginst_d_obs_vars = NULL;
139 /* sampling vector for feeding smarter integrators */
140 static struct Integ_samples_t g_intg_samples = {0L, NULL, NULL};
141
142 /* Integ_system_build locally global var */
143 static struct Integ_system_t *l_isys = NULL;
144 static long l_nstates = 0;
145 static long l_nderivs = 0;
146
147 /* Macros for d. arrays.
148 D_Y_VAR(i) returns the atom that corresponds to the FORTRAN y(i).
149 D_DYDX_VAR(i) returns the atom that corresponds to the FORTRAN ydot(i).
150 D_OBS_VAR(i) returns the atom that corresponds to the ASCEND d.obs[i].
151 */
152 #define D_Y_VAR(i) g_intginst_d_y_vars[(i)-1]
153 #define D_DYDX_VAR(i) g_intginst_d_dydx_vars[(i)-1]
154 #define D_OBS_VAR(i) g_intginst_d_obs_vars[(i)-1]
155
156 /********************************************************************/
157
158 int Asc_IntegCheckStatus(slv_status_t status) {
159 if (status.converged) {
160 return 1;
161 }
162 if (status.diverged) {
163 FPRINTF(stderr, "The derivative system did not converge.\n");
164 FPRINTF(stderr, "Integration will be terminated ");
165 FPRINTF(stderr, "at the end of the current step.\n");
166 return 0;
167 }
168 if (status.inconsistent) {
169 FPRINTF(stderr, "A numerical inconsistency was discovered ");
170 FPRINTF(stderr, "while calculating derivatives.");
171 FPRINTF(stderr, "Integration will be terminated at the end of ");
172 FPRINTF(stderr, "the current step.\n");
173 return 0;
174 }
175 if (status.time_limit_exceeded) {
176 FPRINTF(stderr, "The time limit was ");
177 FPRINTF(stderr, "exceeded while calculating derivatives.\n");
178 FPRINTF(stderr, "Integration will be terminated at ");
179 FPRINTF(stderr, "the end of the current step.\n");
180 return 0;
181 }
182 if (status.iteration_limit_exceeded) {
183 FPRINTF(stderr, "The iteration limit was ");
184 FPRINTF(stderr, "exceeded while calculating derivatives.\n");
185 FPRINTF(stderr, "Integration will be terminated at ");
186 FPRINTF(stderr, "the end of the current step.\n");
187 return 0;
188 }
189 if (status.panic) {
190 FPRINTF(stderr, "The user patience limit was ");
191 FPRINTF(stderr, "exceeded while calculating derivatives.\n");
192 FPRINTF(stderr, "Integration will be terminated at ");
193 FPRINTF(stderr, "the end of the current step.\n");
194 return 0;
195 }
196 return 0;
197 }
198
199 /********************************************************************/
200 static FILE *l_y_file = NULL;
201 static FILE *l_obs_file = NULL;
202 static char *l_y_filename = NULL;
203 static char *l_obs_filename = NULL;
204 static int l_print_option = 1; /* default si */
205 static int l_print_fmt = 0; /* default variable */
206
207 FILE *Asc_IntegOpenYFile(void)
208 {
209 if (l_y_filename==NULL) {
210 return NULL;
211 }
212 l_y_file = fopen(l_y_filename,"a+");
213
214 if (l_y_file==NULL) {
215 FPRINTF(ASCERR,
216 "WARNING: (integrate) Unable to open\n\t%s\nfor state output log.\n",
217 l_y_filename);
218 } else {
219 time_t t;
220 t = time((time_t *)NULL);
221 FPRINTF(l_y_file,"DATASET %s", asctime(localtime(&t)));
222 FFLUSH(l_y_file);
223 }
224 return l_y_file;
225 }
226
227 FILE *Asc_IntegOpenObsFile(void)
228 {
229 if (l_obs_filename==NULL) {
230 return NULL;
231 }
232 l_obs_file = fopen(l_obs_filename,"a+");
233 if (l_obs_file==NULL) {
234 FPRINTF(ASCERR,
235 "WARNING: (integrate) Unable to open\n\t%s\nfor observation log.\n",
236 l_obs_filename);
237 } else {
238 time_t t;
239 t = time((time_t *)NULL);
240 FPRINTF(l_obs_file,"DATASET %s", asctime(localtime(&t)));
241 FFLUSH(l_obs_file);
242 }
243 return l_obs_file;
244 }
245
246 FILE *Asc_IntegGetYFile(void)
247 {
248 return l_y_file;
249 }
250
251 FILE *Asc_IntegGetObsFile(void)
252 {
253 return l_obs_file;
254 }
255
256 void Asc_IntegReleaseYFile(void)
257 {
258 l_y_file = NULL;
259 }
260 void Asc_IntegReleaseObsFile(void)
261 {
262 l_obs_file = NULL;
263 }
264 /********************************************************************/
265 /* string column layout: 26 char columns if l_print_fmt else variable */
266 /* numeric format: leading space plus characters from Unit*Value */
267 /* index format: left padded long int */
268 /* headline format space plus - s */
269 #define BCOLSFMT (l_print_fmt ? "%-26s" : "\t%s")
270 #define BCOLNFMT (l_print_fmt ? " %-25s" : "\t%s")
271 #define BCOLIFMT (l_print_fmt ? " %25ld" : "\t%ld")
272 #define BCOLHFMT (l_print_fmt ? " -------------------------" : "\t---")
273
274 void Asc_IntegPrintYHeader(FILE *fp, struct Integ_system_t *blsys)
275 {
276 long i,len, *yip;
277 char *name;
278 struct Instance *in;
279 int si;
280 if (fp==NULL) {
281 return;
282 }
283 if (blsys==NULL) {
284 FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYHeader: called w/o data\n");
285 return;
286 }
287 if (blsys->n_y == 0) {
288 return;
289 }
290 if (blsys->y == NULL) {
291 FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");
292 return;
293 }
294 len = blsys->n_y;
295 yip = blsys->y_id;
296 si = l_print_option;
297 /* output indep var name */
298 /* output dep var names */
299 FPRINTF(fp,"States: (user index) (name) (units)\n");
300 in = var_instance(blsys->x);
301 FPRINTF(fp,"{indvar}"); /* indep id name */
302 name = WriteInstanceNameString(in,g_solvinst_cur);
303 FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
304 ascfree(name);
305 for (i=0; i< len; i++) {
306 in = var_instance(blsys->y[i]);
307 FPRINTF(fp,"{%ld}",yip[i]); /* user id # */
308 name = WriteInstanceNameString(in,g_solvinst_cur);
309 FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
310 ascfree(name);
311 }
312 FPRINTF(fp,BCOLSFMT,"indvar");
313 for (i=0; i < len; i++) {
314 FPRINTF(fp,BCOLIFMT,yip[i]);
315 }
316 FPRINTF(fp,"\n");
317 for (i=0; i <= len; i++) {
318 FPRINTF(fp,BCOLHFMT);
319 }
320 FPRINTF(fp,"\n");
321 }
322 /********************************************************************/
323 void Asc_IntegPrintObsHeader(FILE *fp, struct Integ_system_t *blsys)
324 {
325 long i,len, *obsip;
326 char *name;
327 struct Instance *in;
328 int si;
329 if (fp==NULL) {
330 return;
331 }
332 if (blsys==NULL) {
333 FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsHeader: called w/o data\n");
334 return;
335 }
336 if (blsys->n_obs == 0) {
337 return;
338 }
339 if (blsys->obs == NULL) {
340 FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");
341 return;
342 }
343 len = blsys->n_obs;
344 obsip = blsys->obs_id;
345 si = l_print_option;
346 FPRINTF(fp,"Observations: (user index) (name) (units)\n");
347 /* output indep var name */
348 /* output obs var names */
349 in = var_instance(blsys->x);
350 FPRINTF(fp,"{indvar}"); /* indep id name */
351 name = WriteInstanceNameString(in,g_solvinst_cur);
352 FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
353 ascfree(name);
354 for (i=0; i< len; i++) {
355 in = var_instance(blsys->obs[i]);
356 FPRINTF(fp,"{%ld}",obsip[i]); /* user id # */
357 name = WriteInstanceNameString(in,g_solvinst_cur);
358 FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
359 ascfree(name);
360 }
361 FPRINTF(fp,BCOLSFMT,"indvar");
362 for (i=0; i < len; i++) {
363 FPRINTF(fp,BCOLIFMT,obsip[i]);
364 }
365 FPRINTF(fp,"\n");
366 for (i=0; i <= len; i++) {
367 FPRINTF(fp,BCOLHFMT);
368 }
369 FPRINTF(fp,"\n");
370 }
371
372 /********************************************************************/
373 void Asc_IntegPrintYLine(FILE *fp, struct Integ_system_t *blsys)
374 {
375 long i,len;
376 struct var_variable **vp;
377 int si;
378 if (fp==NULL) {
379 return;
380 }
381 if (blsys==NULL) {
382 FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYLine: called w/o data\n");
383 return;
384 }
385 if (blsys->n_y == 0) {
386 return;
387 }
388 if (blsys->y == NULL) {
389 FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");
390 return;
391 }
392 vp = blsys->y;
393 len = blsys->n_y;
394 si = l_print_option;
395 FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));
396 for (i=0; i < len; i++) {
397 FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));
398 }
399 FPRINTF(fp,"\n");
400 }
401
402 void Asc_IntegPrintObsLine(FILE *fp, struct Integ_system_t *blsys)
403 {
404 long i,len;
405 struct var_variable **vp;
406 int si;
407 if (fp==NULL) {
408 return;
409 }
410 if (blsys==NULL) {
411 FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsLine: called w/o data\n");
412 return;
413 }
414 if (blsys->n_obs == 0) {
415 return;
416 }
417 if (blsys->obs == NULL) {
418 FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");
419 return;
420 }
421 vp = blsys->obs;
422 len = blsys->n_obs;
423 si = l_print_option;
424 FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));
425 for (i=0; i < len; i++) {
426 FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));
427 }
428 FPRINTF(fp,"\n");
429 }
430
431 int Asc_IntegSetYFileCmd(ClientData cdata, Tcl_Interp *interp,
432 int argc, CONST84 char *argv[])
433 {
434 size_t len;
435
436 (void)cdata; /* stop gcc whine about unused parameter */
437
438 if ( argc != 2 ) {
439 FPRINTF(ASCERR, "integrate_set_y_file: called without filename.\n");
440 Tcl_SetResult(interp,
441 "integrate_set_y_file <filename,""> called without arg.",
442 TCL_STATIC);
443 return TCL_ERROR;
444 }
445 if (l_y_filename != NULL) {
446 ascfree(l_y_filename);
447 }
448 len = strlen(argv[1]);
449 if (len >0 ) {
450 l_y_filename = Asc_MakeInitString((int)len);
451 sprintf(l_y_filename,"%s",argv[1]);
452 } else {
453 l_y_filename = NULL;
454 }
455 return TCL_OK;
456 }
457
458 int Asc_IntegSetObsFileCmd(ClientData cdata, Tcl_Interp *interp,
459 int argc, CONST84 char *argv[])
460 {
461 size_t len;
462
463 (void)cdata; /* stop gcc whine about unused parameter */
464
465 if ( argc != 2 ) {
466 FPRINTF(ASCERR, "integrate_set_obs_file: called without filename.\n");
467 Tcl_SetResult(interp,
468 "integrate_set_obs_file <filename,""> called without arg.",
469 TCL_STATIC);
470 return TCL_ERROR;
471 }
472 if (l_obs_filename != NULL) {
473 ascfree(l_obs_filename);
474 }
475 len = strlen(argv[1]);
476 if (len >0 ) {
477 l_obs_filename = Asc_MakeInitString((int)len);
478 sprintf(l_obs_filename,"%s",argv[1]);
479 } else {
480 l_obs_filename = NULL;
481 }
482 return TCL_OK;
483 }
484
485 int Asc_IntegSetFileUnitsCmd(ClientData cdata, Tcl_Interp *interp,
486 int argc, CONST84 char *argv[])
487 {
488 (void)cdata; /* stop gcc whine about unused parameter */
489
490 if ( argc != 2 ) {
491 FPRINTF(ASCERR, "integrate_logunits: called without printoption.\n");
492 Tcl_SetResult(interp,"integrate_logunits <display,si> called without arg.",
493 TCL_STATIC);
494 return TCL_ERROR;
495 }
496 switch (argv[1][0]) {
497 case 's':
498 l_print_option = 1;
499 break;
500 case 'd':
501 l_print_option = 0;
502 break;
503 default:
504 FPRINTF(ASCERR,"integrate_logunits: called with bogus argument.\n");
505 FPRINTF(ASCERR,"logunits remain set to %s.\n",
506 (l_print_option ? "si":"display"));
507 break;
508 }
509 return TCL_OK;
510 }
511
512 int Asc_IntegSetFileFormatCmd(ClientData cdata, Tcl_Interp *interp,
513 int argc, CONST84 char *argv[])
514 {
515 (void)cdata; /* stop gcc whine about unused parameter */
516
517 if ( argc != 2 ) {
518 FPRINTF(ASCERR, "integrate_logformat called without printoption.\n");
519 Tcl_SetResult(interp,
520 "integrate_logformat <fixed,variable> called without arg.",
521 TCL_STATIC);
522 return TCL_ERROR;
523 }
524 switch (argv[1][0]) {
525 case 'f':
526 l_print_fmt = 1;
527 break;
528 case 'v':
529 l_print_fmt = 0;
530 break;
531 default:
532 FPRINTF(ASCERR,"integrate_logformat: called with bogus argument.\n");
533 FPRINTF(ASCERR,"logformat remains set to %s.\n",
534 (l_print_fmt ? "fixed":"variable"));
535 break;
536 }
537 return TCL_OK;
538 }
539
540
541 /********************************************************************/
542 /* stuff for handling a user defined list of step outside the ascend
543 model */
544
545 int Asc_IntegSetXSamples(long n, double *values, dim_type *d)
546 {
547 /* trash the old stuff, regardless.*/
548 g_intg_samples.ns = 0L;
549 if (g_intg_samples.x != NULL) {
550 ascfree(g_intg_samples.x);
551 g_intg_samples.x = NULL;
552 }
553 if (g_intg_samples.d != NULL) {
554 ascfree(g_intg_samples.d);
555 g_intg_samples.d = NULL;
556 }
557 /* if was a reset, return now */
558 if (n <1 || values==NULL) {
559 return 0;
560 }
561 /* store d given or copy and store WildDimension */
562 if (d != NULL) {
563 g_intg_samples.d = d;
564 #if INTDEBUG
565 FPRINTF(ASCERR,"sample received dimen\n");
566 PrintDimen(ASCERR,g_intg_samples.d);
567 #endif
568 } else {
569 g_intg_samples.d = (dim_type *)ascmalloc(sizeof(dim_type));
570 if (g_intg_samples.d == NULL) {
571 FPRINTF(ASCERR,"ERROR: (Asc_IntegSetXSamples) Insufficient memory.\n");
572 return 1;
573 }
574 CopyDimensions(WildDimension(),g_intg_samples.d);
575 #if INTDEBUG
576 FPRINTF(ASCERR,"copy of wild dimen looks like\n");
577 PrintDimen(ASCERR,g_intg_samples.d);
578 #endif
579 }
580 g_intg_samples.ns = n;
581 g_intg_samples.x = values;
582 return 0;
583 }
584
585 int Asc_IntegGetMaxSteps(struct Integ_system_t *sys)
586 {
587 return sys->maxsteps;
588 }
589
590 double Asc_IntegGetStepMax(struct Integ_system_t *sys)
591 {
592 return sys->stepmax;
593 }
594
595 double Asc_IntegGetStepMin(struct Integ_system_t *sys)
596 {
597 return sys->stepmin;
598 }
599
600 double Asc_IntegGetStepZero(struct Integ_system_t *sys)
601 {
602 return sys->stepzero;
603 }
604
605 long Asc_IntegGetNSamples()
606 {
607 return g_intg_samples.ns;
608 }
609
610 double Asc_IntegGetXSamplei(long i)
611 {
612 if (i > -1 && i < g_intg_samples.ns) {
613 return g_intg_samples.x[i];
614 } else {
615 FPRINTF(ASCERR,
616 "WARNING: (Asc_IntegGetXSamplei) Undefined xsample %ld."
617 " Returning 0.0.\n",
618 i);
619 return (double)0.0;
620 }
621 }
622
623 void Asc_IntegSetXSamplei(long i,double xi)
624 {
625 if (i > -1 && i < g_intg_samples.ns) {
626 g_intg_samples.x[i] = xi;
627 } else {
628 FPRINTF(ASCERR,
629 "WARNING: (Asc_IntegSetXSamplei) Undefined xsample %ld. Ignored.\n",
630 i);
631 }
632 }
633
634 dim_type *Asc_IntegGetXSamplesDim(void)
635 {
636 return g_intg_samples.d;
637 }
638 /********************************************************************/
639 /* callbacks for manipulating a user defined list of steps */
640 int Asc_IntegGetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,
641 int argc, CONST84 char *argv[])
642 {
643 static char sval[40]; /* buffer long enough to hold a double printed */
644 struct Units *du = NULL;
645 dim_type *dp;
646 long i,len;
647 double *uvalues = NULL;
648 char *ustring;
649 double *uv;
650 int trydu=0, prec, stat=0;
651
652 (void)cdata; /* stop gcc whine about unused parameter */
653
654 if (( argc < 1 ) || ( argc > 2 )) {
655 Tcl_SetResult(interp,
656 "integrate_get_samples: expected 0 or 1 args [display]",
657 TCL_STATIC);
658 return TCL_ERROR;
659 }
660 if (argc==2) {
661 trydu = 1;
662 if( argv[1][0] != 'd') {
663 Tcl_SetResult(interp, "integrate_get_samples: expected display but got ",
664 TCL_STATIC);
665 Tcl_AppendResult(interp,argv[1],".",SNULL);
666 return TCL_ERROR;
667 }
668 }
669 len = g_intg_samples.ns;
670 dp = g_intg_samples.d;
671 if (len <1) {
672 Tcl_SetResult(interp, "{} {}", TCL_STATIC);
673 return TCL_OK;
674 }
675
676 if (trydu) {
677 uvalues = (double *)ascmalloc(sizeof(double)*len);
678 if (uvalues == NULL) {
679 Tcl_SetResult(interp, "integrate_get_samples: Insufficient memory.",
680 TCL_STATIC);
681 return TCL_ERROR;
682 }
683 ustring = Asc_UnitDimString(dp,0); /* get display unit string */
684 du = (struct Units *)LookupUnits(ustring);
685 if (du == NULL) {
686 /* a very bizarre thing to happen,
687 * since Asc_UnitDimString just made it
688 */
689 stat = 1;
690 } else {
691 stat = 0;
692 uv = uvalues;
693 for (i=0; i < len; i++) {
694 stat = Asc_UnitConvert(du,g_intg_samples.x[i],uv,1);
695 if (stat) {
696 break;
697 }
698 uv++;
699 }
700 }
701 if (stat) {
702 ascfree(uvalues);
703 }
704 }
705 /* if display not wanted or display -> conversion error */
706 if (!trydu || stat) {
707 uvalues = g_intg_samples.x;
708 ustring = Asc_UnitDimString(dp,1); /* get si unit string */
709 }
710
711 Tcl_AppendResult(interp,"{",ustring,"} {",SNULL);
712 prec = Asc_UnitGetCPrec();
713 len--; /* last one is a special case */
714 for (i=0; i<len; i++) {
715 /* print number and a blank */
716 sprintf(sval,"%.*g ",prec,uvalues[i]);
717 Tcl_AppendResult(interp,sval,SNULL);
718 }
719 sprintf(sval,"%.*g",prec,uvalues[len]);
720 Tcl_AppendResult(interp,sval,"}",SNULL);
721 if (trydu && !stat) {
722 ascfree(uvalues);
723 }
724 return TCL_OK;
725 }
726
727 int Asc_IntegSetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,
728 int argc, CONST84 char *argv[])
729 {
730 struct Units *du = NULL;
731 dim_type *dp;
732 dim_type *mydp;
733 long i,len;
734 double *uvalues = NULL;
735 double *uv;
736 int stat;
737
738 (void)cdata; /* stop gcc whine about unused parameter */
739
740 if (argc == 1) {
741 Asc_IntegSetXSamples(0L,NULL,NULL);
742 return TCL_OK;
743 }
744
745 if (argc <4) {
746 Tcl_SetResult(interp,
747 "Syntax: integrate_set_samples"
748 " <units> <value [value...] value>",
749 TCL_STATIC);
750 FPRINTF(ASCERR,"ERROR: integrate_set_samples needs at least 3 args.");
751 return TCL_ERROR;
752 }
753
754 du = (struct Units *)LookupUnits(argv[1]);
755 if (du == NULL) {
756 Tcl_SetResult(interp, "integrate_set_samples: first arg not valid units.",
757 TCL_STATIC);
758 return TCL_ERROR;
759 }
760 dp = (dim_type *)UnitsDimensions(du);
761 #if INTDEBUG
762 FPRINTF(ASCERR,"user dimen looks like\n");
763 PrintDimen(ASCERR,dp);
764 #endif
765 mydp = (dim_type *)ascmalloc(sizeof(dim_type)); /* we are giving away */
766 if (mydp == NULL) {
767 Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",
768 TCL_STATIC);
769 return TCL_ERROR;
770 }
771 CopyDimensions(dp,mydp);
772 #if INTDEBUG
773 FPRINTF(ASCERR,"copy of user dimen looks like\n");
774 PrintDimen(ASCERR,mydp);
775 #endif
776
777 len = argc -2;
778 uvalues = (double *)ascmalloc(sizeof(double) * len); /* we are giving away */
779 if (uvalues==NULL) {
780 Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",
781 TCL_STATIC);
782 ascfree(mydp);
783 return TCL_ERROR;
784 }
785 stat = 0;
786 uv = uvalues;
787 for (i=0; i<len; i++) {
788 if(Tcl_GetDouble(interp,argv[i+2],uv)!=TCL_OK) {
789 stat = 1;
790 break;
791 }
792 stat = Asc_UnitConvert(du,*uv,uv,0);
793 if (stat) {
794 break;
795 }
796 uv++;
797 }
798 Tcl_ResetResult(interp);
799 if (stat) {
800 Tcl_SetResult(interp, "integrate_set_samples: Invalid value given. (",
801 TCL_STATIC);
802 Tcl_AppendResult(interp,argv[i+2],")",SNULL);
803 ascfree(uvalues);
804 ascfree(mydp);
805 return TCL_ERROR;
806 }
807 if (Asc_IntegSetXSamples(len,uvalues,mydp)) {
808 Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory.",
809 TCL_STATIC);
810 ascfree(uvalues);
811 ascfree(mydp);
812 return TCL_ERROR;
813 }
814 return TCL_OK;
815 }
816
817 /********************************************************************/
818
819 int Asc_IntegInstIntegrable(struct Instance *inst,
820 enum Integrator_type integrator)
821 {
822 switch (integrator) {
823 case BLSODE :
824 if (inst == NULL) {
825 return 0;
826 }
827 return 1;
828 case UNKNOWN:
829 FPRINTF(stderr, "UNKNOWN integrator is not supported.\n");
830 return 0;
831 default :
832 FPRINTF(stderr, "ERRONEOUS integrator is not supported.\n");
833 return 0;
834 }
835 }
836
837 int Asc_IntegInstIntegrableCmd(ClientData cdata,Tcl_Interp *interp,
838 int argc, CONST84 char *argv[])
839 {
840 struct Instance *i=NULL;
841 enum Integrator_type integrator=UNKNOWN;
842 int result=0; /* 0 = FALSE; 1 = TRUE */
843
844 (void)cdata; /* stop gcc whine about unused parameter */
845
846 if ( argc != 3 ) {
847 Tcl_SetResult(interp, "integrate_able <solver,current,search> <lsode>",
848 TCL_STATIC);
849 return TCL_ERROR;
850 }
851
852 if (strncmp(argv[1],"solver",3)==0) {
853 i = g_solvinst_cur;
854 } else {
855 if (strncmp(argv[1],"search",3)==0) {
856 i = g_search_inst;
857 } else {
858 if (strncmp(argv[1],"current",3)==0) {
859 i = g_curinst;
860 } else {
861 Tcl_SetResult(interp,
862 "integrate_able: arg 1 is current, search, or solver",
863 TCL_STATIC);
864 return TCL_ERROR;
865 }
866 }
867 }
868
869 if (!i) {
870 Tcl_SetResult(interp, "0", TCL_STATIC);
871 FPRINTF(ASCERR,"NULL instance sent to integrate_able.\n");
872 return TCL_OK;
873 }
874
875 if (strncmp(argv[2],"blsode",3)==0) {
876 integrator = BLSODE;
877 }
878 /* someday will have an ivp homegrown. in anycase, ivp dis/en-ables btn */
879 if (strncmp(argv[2],"ivp",3)==0) {
880 integrator = UNKNOWN;
881 }
882
883 result = Asc_IntegInstIntegrable(i, integrator);
884 if (result) {
885 Tcl_SetResult(interp, "1", TCL_STATIC);
886 } else {
887 Tcl_SetResult(interp, "0", TCL_STATIC);
888 }
889 return TCL_OK;
890 }
891
892 /*** derivative parts **** x *****************************************/
893
894 double Asc_IntegGetDX() {
895 return var_value(g_intginst_d_x);
896 }
897
898 void Asc_IntegSetDX( double value) {
899 var_set_value(g_intginst_d_x, value);
900 print_debug("set_d_x = %g\n", value);
901 }
902
903 /*** derivative parts **** y *****************************************/
904
905 double *Asc_IntegGetDY(double *y) {
906 long i;
907
908 if (y==NULL) {
909 y = (double *)asccalloc((g_intginst_d_n_eq+1), sizeof(double));
910 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
911 }
912
913 for (i=1; i<=g_intginst_d_n_eq; i++) {
914 y[i-1] = var_value(D_Y_VAR(i));
915 print_debug("*get_d_y[%ld] = ", i);
916 print_debug("%g\n", y[i-1]);
917 }
918 return y;
919 }
920
921 void Asc_IntegSetDY(double *y) {
922 long i;
923
924 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
925
926 for (i=1; i<=g_intginst_d_n_eq; i++) {
927 var_set_value(D_Y_VAR(i),y[i-1]);
928 print_debug("*set_d_y[%ld] = ", i);
929 print_debug("%g\n", y[i-1]);
930 }
931 }
932
933 /*** derivative parts **** dydx *****************************************/
934
935 double *Asc_IntegGetDDydx(double *dydx) {
936 long i;
937
938 if (dydx==NULL) {
939 dydx = (double *)asccalloc((g_intginst_d_n_eq+1), sizeof(double));
940 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
941 }
942
943 for (i=1; i<=g_intginst_d_n_eq; i++) {
944 dydx[i-1] = var_value(D_DYDX_VAR(i));
945 print_debug("*get_d_dydx[%ld] = ", i);
946 print_debug("%g\n", dydx[i-1]);
947 }
948 return dydx;
949 }
950
951 void Asc_IntegSetDDydx(double *dydx) {
952 long i;
953
954 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
955
956 for (i=1; i<=g_intginst_d_n_eq; i++) {
957 var_set_value(D_DYDX_VAR(i),dydx[i-1]);
958 print_debug("*set_d_dydx[%ld] = ", i);
959 print_debug("%g\n", dydx[i-1]);
960 }
961 }
962
963 /**** derivative parts * d.obs ******************************************/
964
965 /*
966 This function takes the inst in the solver and returns the vector of
967 observation variables that are located in the submodel d.obs array.
968 */
969 double *Asc_IntegGetDObs(double *obsi) {
970 long i;
971
972 if (obsi==NULL) {
973 obsi = (double *)asccalloc((g_intginst_d_n_obs+1),sizeof(double));
974 }
975
976 /* C obsi[0] <==> ascend d.obs[1] */
977
978 for (i=1; i<=g_intginst_d_n_obs; i++) {
979 obsi[i-1] = var_value(D_OBS_VAR(i));
980 print_debug("*get_d_obs[%ld] = ", i);
981 print_debug("%g\n", obsi[i-1]);
982 }
983 return obsi;
984 }
985
986 /* All the KEEPIP stuff used to be here. its gone to cvs land */
987
988 /*
989 * takes the type of integrator and start and finish index and calls the
990 * appropriate integrator
991 */
992 static void Integ_Solve( enum Integrator_type integrator,
993 int start_index, int finish_index,
994 struct Integ_system_t *blsys) {
995 switch (integrator) {
996 case BLSODE:
997 Asc_BLsodeIntegrate(g_intgsys_cur,start_index, finish_index, blsys);
998 return;
999 default:
1000 FPRINTF(stderr, "The requested integrator is ");
1001 FPRINTF(stderr, "not currently available\n.");
1002 return;
1003 }
1004 }
1005
1006 static double **MakeDenseMatrix(int nrows, int ncols)
1007 {
1008 int c;
1009 double **result;
1010 assert(nrows>0);
1011 assert(ncols>0);
1012 result = (double **)ascmalloc(nrows*sizeof(double *));
1013 for (c=0;c<nrows;c++) {
1014 result[c] = (double *)asccalloc(ncols,sizeof(double));
1015 }
1016 return result;
1017 }
1018
1019 static void DestroyDenseMatrix(double **matrix,int nrows)
1020 {
1021 int c;
1022 if (matrix) {
1023 for (c=0;c<nrows;c++) {
1024 if (matrix[c]) {
1025 ascfree((char *)matrix[c]);
1026 }
1027 }
1028 ascfree((char *)matrix);
1029 }
1030 }
1031
1032 /*
1033 * needs work. Assumes struct Instance* and struct var_variable*
1034 * are synonymous, which demonstrates the need for a method to take
1035 * an instance and ask the solvers for its global or local index
1036 * if var and inst are decoupled.
1037 */
1038 static
1039 int Integ_SetUpDiffs_BLsode(struct Integ_system_t *blsys) {
1040 long n_eqns;
1041 unsigned long nch,i;
1042
1043 struct var_variable **vp;
1044 int *ip;
1045
1046 g_intg_diff.n_eqns = n_eqns = blsys->n_y;
1047 g_intg_diff.input_indices = (int *)asccalloc(n_eqns, sizeof(int));
1048 g_intg_diff.output_indices = (int *)asccalloc(n_eqns, sizeof(int));
1049 g_intg_diff.y_vars = NULL;
1050 g_intg_diff.ydot_vars = NULL;
1051 g_intg_diff.dydx_dx = MakeDenseMatrix(n_eqns,n_eqns);
1052
1053
1054 /*
1055 * Let us now process what we consider *inputs* to the problem as
1056 * far as ASCEND is concerned; i.e. the state vars or the y_vars's
1057 * if you prefer.
1058 */
1059 nch = n_eqns;
1060 vp = g_intg_diff.y_vars =
1061 (struct var_variable **)ascmalloc((nch+1)*sizeof(struct var_variable *));
1062 ip = g_intg_diff.input_indices;
1063 for (i=0;i<nch;i++) {
1064 *vp = (struct var_variable *)blsys->y[i];
1065 *ip = var_sindex(*vp);
1066 vp++;
1067 ip++;
1068 }
1069 *vp = NULL; /* terminate */
1070
1071 /*
1072 * Let us now go for the outputs, ie the derivative terms.
1073 */
1074 vp = g_intg_diff.ydot_vars =
1075 (struct var_variable **)ascmalloc((nch+1)*sizeof(struct var_variable *));
1076 ip = g_intg_diff.output_indices;
1077 for (i=0;i<nch;i++) {
1078 *vp = (struct var_variable *)blsys->ydot[i];
1079 *ip = var_sindex(*vp);
1080 vp++; /* dont assume that a var is synonymous with */
1081 ip++; /* an Instance; that might/will change soon */
1082 }
1083 *vp = NULL; /* terminate */
1084
1085 return 0;
1086 }
1087
1088 /* Build an analytic jacobian for solving the state system */
1089 /*
1090 * This necessarily ugly piece of code attempts to create a unique
1091 * list of relations that explicitly contain the variables in the
1092 * given input list. The utility of this information is that we know
1093 * exactly which relations must be differentiated, to fill in the
1094 * df/dy matrix. If the problem has very few derivative terms, this will
1095 * be of great savings. If the problem arose from the discretization of
1096 * a pde, then this will be not so useful. The decision wether to use
1097 * this function or to simply differentiate the entire relations list
1098 * must be done before calling this function. Final Note: the callee
1099 * owns the array, but not the array elements.
1100 */
1101
1102 #define AVG_NUM_INCIDENT 4
1103
1104 /* Returns STATEFLAG child value if we have a correct dynamic variable.
1105 If we do not, returns 0. Doesn't check solver_var status.
1106 if index pointer is given, it is stuffed too.
1107 if return is 0, ignore index.
1108 to be correct either type and index are assigned integers or
1109 type is -1. */
1110 static long DynamicVarInfo(struct var_variable *v,long *index)
1111 {
1112 struct Instance *c, *d, *i;
1113 i = var_instance(v);
1114 assert(i!=NULL);
1115 c = ChildByChar(i,STATEFLAG);
1116 d = ChildByChar(i,STATEINDEX);
1117 /* lazy evaluation is important in the following if */
1118 if( c == NULL ||
1119 d == NULL ||
1120 InstanceKind(c) != INTEGER_INST ||
1121 InstanceKind(d) != INTEGER_INST ||
1122 !AtomAssigned(c) ||
1123 (!AtomAssigned(d) && GetIntegerAtomValue(c) != -1L)
1124 ) {
1125 return 0L;
1126 }
1127 if (index != NULL) {
1128 *index = GetIntegerAtomValue(d);
1129 }
1130 return GetIntegerAtomValue(c);
1131 }
1132
1133 /* returns the pointer if we have a correct observation variable.
1134 If long is passed in, long will have the index value.
1135 Vars with UNDEFINED observation flags are 'incorrect.'
1136 */
1137 static struct var_variable *ObservationVar(struct var_variable *v, long *index)
1138 {
1139 struct Instance *c,*i;
1140 i = var_instance(v);
1141 assert(i!=NULL);
1142 c = ChildByChar(i,OBSINDEX);
1143 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1144 return NULL;
1145 }
1146 if (index != NULL) {
1147 *index = GetIntegerAtomValue(c);
1148 }
1149 return v;
1150 }
1151
1152 /* take the variable instance and set its obsid to index.
1153 id must be already assigned at least once. */
1154 static void Integ_SetObsId(struct var_variable *v, long index)
1155 {
1156 struct Instance *c, *i;
1157 i = var_instance(v);
1158 assert(i!=NULL);
1159 c = ChildByChar(i,OBSINDEX);
1160 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1161 return;
1162 }
1163 SetIntegerAtomValue(c,index,0);
1164 }
1165
1166 /*
1167 * Sorts the interesting vars into our wonderful little lists.
1168 * Dislikes null input intensely.
1169 */
1170 static void Integ_classify_vars(struct var_variable *var)
1171 {
1172 struct Integ_var_t *info;
1173 long type,index;
1174 var_filter_t vfilt;
1175
1176 assert(var != NULL && var_instance(var)!=NULL );
1177 vfilt.matchbits = VAR_SVAR;
1178 vfilt.matchvalue = VAR_SVAR;
1179 if( var_apply_filter(var,&vfilt) ) {
1180 type = DynamicVarInfo(var,&index);
1181 if ( type != 0) {
1182 #if INTEG_DEBUG
1183 var_write_name(g_solvsys_cur,var,ASCERR);
1184 FPRINTF(ASCERR," type = %ld\n",type);
1185 #endif
1186 /* ignore algebraics type == 0 */
1187 info = (struct Integ_var_t *)ascmalloc(sizeof(struct Integ_var_t));
1188 if (info == NULL) {
1189 FPRINTF(ASCERR,
1190 "ERROR: (Integ_system_build) Insufficient memory.\n");
1191 return;
1192 }
1193 info->type = type;
1194 info->index = index;
1195 info->i = var;
1196 if (type > 0L) {
1197 /* state or derivative */
1198 gl_append_ptr(l_isys->dynvars,(POINTER)info);
1199 if (type == 1) {
1200 l_nstates++;
1201 }
1202 if (type == 2) {
1203 l_nderivs++;
1204 }
1205 } else {
1206 if (type == -1L) {
1207 /* independent */
1208 gl_append_ptr(l_isys->indepvars,(POINTER)info);
1209 } else {
1210 /* probably should whine about unknown type here */
1211 ascfree(info);
1212 }
1213 }
1214 info = NULL;
1215 }
1216 if ( ObservationVar(var,&index) != NULL && index > 0L) {
1217 info = (struct Integ_var_t *)ascmalloc(sizeof(struct Integ_var_t));
1218 if (info == NULL) {
1219 FPRINTF(ASCERR,
1220 "ERROR: (Integ_system_build) Insufficient memory.\n");
1221 return;
1222 }
1223 info->type = 0L;
1224 info->index = index;
1225 info->i = var;
1226 gl_append_ptr(l_isys->obslist,(POINTER)info);
1227 info = NULL;
1228 }
1229 }
1230 }
1231
1232
1233 /* compares observation structs. NULLs should end up at far end. */
1234 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2)
1235 {
1236 if (v1 == NULL) {
1237 return 1;
1238 }
1239 if (v2 == NULL) {
1240 return -1;
1241 }
1242 if (v1->index > v2->index) {
1243 return 1;
1244 }
1245 if (v1->index == v2->index) {
1246 return 0;
1247 }
1248 return -1;
1249 }
1250
1251 /*
1252 Compares dynamic vars structs. NULLs should end up at far end.
1253 List should be sorted primarily by index and then by type, in order
1254 of increasing value of both.
1255 */
1256 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2)
1257 {
1258 if (v1 == NULL) {
1259 return 1;
1260 }
1261 if (v2 == NULL) {
1262 return -1;
1263 }
1264 if (v1->index > v2->index) {
1265 return 1;
1266 }
1267 if (v1->index != v2->index) {
1268 return -1;
1269 }
1270 if (v1->type > v2->type) {
1271 return 1;
1272 }
1273 return -1;
1274 }
1275
1276 /* trash nonnull contents of a system_t . does ree the pointer sent.*/
1277 static void Integ_system_destroy(struct Integ_system_t *sys)
1278 {
1279 if (sys==NULL) {
1280 return;
1281 }
1282 if (sys->states != NULL) {
1283 gl_destroy(sys->states);
1284 }
1285 if (sys->derivs != NULL) {
1286 gl_destroy(sys->derivs);
1287 }
1288 if (sys->dynvars != NULL) {
1289 gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
1290 }
1291 if (sys->obslist != NULL) {
1292 gl_free_and_destroy(sys->obslist); /* and obslist */
1293 }
1294 if (sys->indepvars != NULL) {
1295 gl_free_and_destroy(sys->indepvars); /* and indepvars */
1296 }
1297 if (sys->y != NULL && !sys->ycount) {
1298 ascfree(sys->y);
1299 }
1300 if (sys->ydot != NULL && !sys->ydotcount) {
1301 ascfree(sys->ydot);
1302 }
1303 if (sys->y_id != NULL) {
1304 ascfree(sys->y_id);
1305 }
1306 if (sys->obs != NULL && !sys->obscount) {
1307 ascfree(sys->obs);
1308 }
1309 if (sys->obs_id != NULL) {
1310 ascfree(sys->obs_id);
1311 }
1312 ascfree(sys);
1313 }
1314
1315 static
1316 void IntegInitSymbols(void)
1317 {
1318 STATEFLAG = AddSymbol("ode_type");
1319 STATEINDEX = AddSymbol("ode_id");
1320 OBSINDEX = AddSymbol("obs_id");
1321 }
1322
1323 /*
1324 * Returns a pointer to an ok struct Integ_system_t. Returns NULL
1325 * if one cannot be built from the instance given.
1326 * Diagnostics will be printed.
1327 * The pointer returned will also be stored in l_isys if someone
1328 * in Integrators.c needs it. Anyone who takes a copy of l_isys should
1329 * then set l_isys to NULL and keep track of the sys from there.
1330 * On exit the gl_list_t * in the returned pointer will all be NULL.
1331 */
1332 static struct Integ_system_t *Integ_system_build(CONST slv_system_t slvsys)
1333 {
1334 struct Integ_system_t *sys;
1335 struct Integ_var_t *v1,*v2;
1336 struct var_variable **vlist;
1337 long half,i,len,vlen;
1338 int happy=1;
1339
1340 if (slvsys==NULL) {
1341 return NULL;
1342 }
1343 sys=(struct Integ_system_t *)asccalloc(1,sizeof(struct Integ_system_t));
1344 if (sys==NULL) {
1345 return sys;
1346 }
1347 l_isys = sys;
1348 IntegInitSymbols();
1349
1350 /* collect potential states and derivatives */
1351 sys->indepvars = gl_create(10L); /* t var info */
1352 sys->dynvars = gl_create(200L); /* y ydot var info */
1353 sys->obslist = gl_create(100L); /* obs info */
1354 if (sys->dynvars == NULL ||
1355 sys->obslist == NULL ||
1356 sys->indepvars == NULL ) {
1357 FPRINTF(ASCERR,"ERROR: (Integ_system_build) Insufficient memory.\n");
1358 Integ_system_destroy(sys);
1359 l_isys = NULL;
1360 return l_isys;
1361 }
1362 l_nstates = l_nderivs = 0;
1363 /* visit all the slv_system_t master var lists to collect vars */
1364 /* find the vars mostly in this one */
1365 vlist = slv_get_master_var_list(slvsys);
1366 vlen = slv_get_num_master_vars(slvsys);
1367 for (i=0;i<vlen;i++) {
1368 Integ_classify_vars(vlist[i]);
1369 }
1370 /* probably nothing here, but gotta check. */
1371 vlist = slv_get_master_par_list(slvsys);
1372 vlen = slv_get_num_master_pars(slvsys);
1373 for (i=0;i<vlen;i++) {
1374 Integ_classify_vars(vlist[i]);
1375 }
1376 /* might find t here */
1377 vlist = slv_get_master_unattached_list(slvsys);
1378 vlen = slv_get_num_master_unattached(slvsys);
1379 for (i=0;i<vlen;i++) {
1380 Integ_classify_vars(vlist[i]);
1381 }
1382
1383 /* check the sanity of the independent variable */
1384 len = gl_length(sys->indepvars);
1385 if (!len) {
1386 FPRINTF(ASCERR,
1387 "ERROR: (Integ_system_build) No independent variable found.\n");
1388 Integ_system_destroy(sys);
1389 l_isys = NULL;
1390 return l_isys;
1391 }
1392 if (len > 1) {
1393 char *name;
1394 FPRINTF(ASCERR,
1395 "ERROR: (Integ_system_build) Excess %ld independent variables found:\n",
1396 len);
1397 for (i=1; i <=len;i++) {
1398 v1 = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
1399 if (v1 != NULL) {
1400 name = WriteInstanceNameString(var_instance(v1->i),g_solvinst_cur);
1401 if (name !=NULL) {
1402 FPRINTF(ASCERR,"\t\n");
1403 ascfree(name);
1404 name = NULL;
1405 }
1406 }
1407 }
1408 FPRINTF(ASCERR, "Set the %s on all but one of these to %s >= 0.\n",
1409 SCP(STATEFLAG),SCP(STATEFLAG));
1410 Integ_system_destroy(sys);
1411 l_isys = NULL;
1412 return l_isys;
1413 } else {
1414 v1 = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
1415 sys->x = v1->i;
1416 }
1417 /* check sanity of state and var lists */
1418 len = gl_length(sys->dynvars);
1419 half = len/2;
1420 if (len % 2 || len == 0L || l_nstates != l_nderivs ) {
1421 /* list length must be even for vars to pair off */
1422 FPRINTF(ASCERR,"ERROR: (Integ_system_build) n_y != n_ydot\n");
1423 FPRINTF(ASCERR,"\tor no dynamic vars found. Fix your indexing.\n");
1424 Integ_system_destroy(sys);
1425 l_isys = NULL;
1426 return l_isys;
1427 }
1428 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
1429 if (gl_fetch(sys->dynvars,len)==NULL) {
1430 FPRINTF(ASCERR,"ERROR: (Integ_system_build) Mysterious NULL found.\n");
1431 FPRINTF(ASCERR,"\tPlease Report this to %s.\n",ASC_BIG_BUGMAIL);
1432 Integ_system_destroy(sys);
1433 l_isys = NULL;
1434 return l_isys;
1435 }
1436 sys->states = gl_create(half); /* state vars Integ_var_t references */
1437 sys->derivs = gl_create(half); /* derivative var atoms */
1438 for (i=1;i < len; i+=2) {
1439 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
1440 v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
1441 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
1442 FPRINTF(ASCERR,"ERROR: (Integ_system_build) Mistyped or misindexed\n");
1443 FPRINTF(ASCERR,
1444 "\tdynamic variables: (%s = %ld,%s = %ld),(%s = %ld,%s = %ld).\n",
1445 SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
1446 SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index);
1447 happy=0;
1448 break;
1449 } else {
1450 gl_append_ptr(sys->states,(POINTER)v1);
1451 gl_append_ptr(sys->derivs,(POINTER)v2->i);
1452 }
1453 }
1454 if (!happy) {
1455 Integ_system_destroy(sys);
1456 l_isys = NULL;
1457 return l_isys;
1458 }
1459 sys->n_y = half;
1460 sys->y = (struct var_variable **)
1461 ascmalloc(sizeof(struct var_variable *)*half);
1462 sys->y_id = (long *)ascmalloc(sizeof(long)*half);
1463 sys->ydot = (struct var_variable **)
1464 ascmalloc(sizeof(struct var_variable *)*half);
1465 if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
1466 FPRINTF(ASCERR,"ERROR: (Integ_system_build) Insufficient memory.\n");
1467 Integ_system_destroy(sys);
1468 l_isys = NULL;
1469 return l_isys;
1470 }
1471 for (i = 1; i <= half; i++) {
1472 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
1473 sys->y[i-1] = v1->i;
1474 sys->y_id[i-1] = v1->index;
1475 sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
1476 }
1477
1478 /* reindex observations. sort if the user mostly numbered. take
1479 * natural order if user just booleaned.
1480 */
1481 len = gl_length(sys->obslist);
1482 /* we shouldn't be seeing NULL here ever except if malloc fail. */
1483 if (len > 1L) {
1484 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
1485 /* half != 0 now because we didn't collect 0 indexed vars */
1486 for (i=2; i <= len; i++) {
1487 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
1488 /* change seen. sort and go on */
1489 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
1490 break;
1491 }
1492 }
1493 }
1494 for (i = half = 1; i <= len; i++) {
1495 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
1496 if (v2==NULL) {
1497 /* we shouldn't be seeing NULL here ever except if malloc fail. */
1498 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
1499 } else {
1500 Integ_SetObsId(v2->i,half);
1501 v2->index = half++;
1502 }
1503 }
1504
1505 /* obslist now uniquely indexed, no nulls */
1506 /* make into arrays */
1507 half = gl_length(sys->obslist);
1508 sys->obs = (struct var_variable **)
1509 ascmalloc(sizeof(struct var_variable *)*half);
1510 sys->obs_id = (long *)ascmalloc(sizeof(long)*half);
1511 if ( sys->obs==NULL || sys->obs_id==NULL) {
1512 FPRINTF(ASCERR,"ERROR: (Integ_system_build) Insufficient memory.\n");
1513 Integ_system_destroy(sys);
1514 l_isys = NULL;
1515 return l_isys;
1516 }
1517 sys->n_obs = half;
1518 for (i = 1; i <= half; i++) {
1519 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
1520 sys->obs[i-1] = v2->i;
1521 sys->obs_id[i-1] = v2->index;
1522 }
1523
1524 /* don't need the gl_lists now that we have arrays for everyone */
1525 gl_destroy(sys->states);
1526 gl_destroy(sys->derivs);
1527 gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
1528 gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
1529 gl_free_and_destroy(sys->obslist); /* and obslist */
1530 sys->states = NULL;
1531 sys->derivs = NULL;
1532 sys->indepvars = NULL;
1533 sys->dynvars = NULL;
1534 sys->obslist = NULL;
1535 return l_isys;
1536 }
1537
1538 /* boolean return. 0 = bad 1 = ok. BLSODE required setup stuff. */
1539 static int Integ_setup_blsode(enum Integrator_type integrator)
1540 {
1541 struct Integ_system_t *blsys=NULL;
1542
1543 g_intgsys_cur = g_solvsys_cur;
1544 if (g_intgsys_cur == NULL) {
1545 FPRINTF(stderr, "g_intgsys_cur not correctly assigned.\n");
1546 return 0;
1547 }
1548 /* verify integrator type ok. always passes for nonNULL inst. */
1549 if (!(Asc_IntegInstIntegrable(g_solvinst_cur, integrator))) {
1550 FPRINTF(stderr, "System is of wrong type to be ");
1551 FPRINTF(stderr, "integrated using BLSODE.\n");
1552 return 0;
1553 }
1554
1555 /* this is a lie, but we will keep it.
1556 We handle any linsol/linsolqr based solver. */
1557 if (strcmp(slv_solver_name(slv_get_selected_solver(g_intgsys_cur)),"QRSlv")
1558 != 0) {
1559 FPRINTF(stderr, "QRSlv must be selected before integration.\n");
1560 return 0;
1561 }
1562
1563 blsys = Integ_system_build(g_solvsys_cur);
1564 if (blsys == NULL) {
1565 FPRINTF(ASCERR,"Unable to build blsode system.\n");
1566 return 0;
1567 }
1568
1569 g_intginst_d_n_eq = blsys->n_y;
1570 g_intginst_d_n_obs = blsys->n_obs;
1571
1572 g_intginst_d_x = blsys->x;
1573 g_intginst_d_y_vars = blsys->y;
1574 g_intginst_d_dydx_vars = blsys->ydot;
1575 g_intginst_d_obs_vars = blsys->obs;
1576 blsys->ycount++;
1577 blsys->ydotcount++;
1578 blsys->obscount++;
1579
1580 print_debugstring("After values\n");
1581 return 1;
1582 }
1583
1584 /*
1585 * takes the type of integrator and sets up the global variables into
1586 * the current integration instance.
1587 */
1588 static int Integ_setup(enum Integrator_type integrator,
1589 long i0, long i1,
1590 double dt0,double dtmin,double dtmax,
1591 int moststeps)
1592 {
1593 long nstep;
1594 unsigned long start_index=0, finish_index=0;
1595 struct Integ_system_t *blsys;
1596
1597 switch (integrator) {
1598 case BLSODE:
1599 if (!Integ_setup_blsode(integrator)) {
1600 return 0;
1601 }
1602 blsys = l_isys;
1603 l_isys = NULL;
1604 break;
1605 default:
1606 FPRINTF(stderr, "The requested type of integrator is ");
1607 FPRINTF(stderr, "not currently available\n.");
1608 return 0;
1609 } /* END of integrator CASE */
1610
1611 switch (integrator) {
1612 case BLSODE:
1613 nstep = Asc_IntegGetNSamples()-1;
1614 /* check for at least 2 steps and dimensionality of x vs steps here */
1615 break;
1616 default:
1617 nstep = 0;
1618 break;
1619 }
1620 if (i0<0 || i1 <0) {
1621 FPRINTF(stdout, "An integration interval had been defined ");
1622 FPRINTF(stdout, "from x[0] to x[%li].\n", nstep);
1623 FPRINTF(stdout, "Enter start index: ");
1624 start_index = (int)readlong((int)start_index);
1625
1626 if (start_index >= (unsigned long)nstep) {
1627 FPRINTF(stderr, "ERROR: Start point < 0 OR start point >= %li.\n",nstep);
1628 Integ_system_destroy(blsys);
1629 return 0;
1630 }
1631
1632 FPRINTF(stdout, "Enter finish index: ");
1633 finish_index = readlong((int)finish_index);
1634
1635 if (finish_index > (unsigned long)nstep) {
1636 FPRINTF(stderr, "ERROR: finish point < 0 OR finish point > %li.\n",nstep);
1637 FPRINTF(stderr, " finish point = %lu\n",finish_index);
1638 Integ_system_destroy(blsys);
1639 return 0;
1640 }
1641 } else {
1642 start_index=i0;
1643 finish_index =i1;
1644 if (start_index >= (unsigned long)nstep) {
1645 FPRINTF(stderr, "ERROR: Start point < 0 OR start point >= %li.\n",nstep);
1646 FPRINTF(stderr, " Start point = %lu\n",start_index);
1647 Integ_system_destroy(blsys);
1648 return 0;
1649 }
1650 if (finish_index > (unsigned long)nstep) {
1651 FPRINTF(stderr, "ERROR: finish point < 0 OR finish point > %li.\n",nstep);
1652 FPRINTF(stderr, " finish point = %lu\n",finish_index);
1653 Integ_system_destroy(blsys);
1654 return 0;
1655 }
1656 }
1657 if ((finish_index <= start_index) || (start_index >= finish_index)) {
1658 FPRINTF(stderr, "ERROR: Finish point <= start point OR ");
1659 FPRINTF(stderr, "start point >= finish point.\n");
1660 FPRINTF(stderr, " Start point = %lu\n",start_index);
1661 FPRINTF(stderr, " finish point = %lu\n",finish_index);
1662 Integ_system_destroy(blsys);
1663 return 0;
1664 }
1665 blsys->maxsteps = moststeps;
1666 blsys->stepzero = dt0;
1667 blsys->stepmin = dtmin;
1668 blsys->stepmax = dtmax;
1669
1670 switch (integrator) {
1671 case BLSODE:
1672 /*
1673 * Set up the Jacobian for the system. Does not check
1674 * This assumes a clean sytem. Shut down will ensure a
1675 * clean system at the end.
1676 */
1677 if (Integ_SetUpDiffs_BLsode(blsys)) {
1678 return 0;
1679 }
1680 break;
1681 default:
1682 break;
1683 }
1684 Integ_Solve(integrator, start_index, finish_index,blsys);
1685 Integ_system_destroy(blsys);
1686 return 0;
1687 }
1688
1689 /*
1690 * switches:
1691 * -engine $name
1692 * -i0 $stepindex
1693 * -i1 $stepindex
1694 * -dt0 $initstepsize
1695 * -dtmin $minstep
1696 * -dtmax $maxstep
1697 * -moststeps $moststeps
1698 */
1699 int Asc_IntegSetupCmd(ClientData cdata,Tcl_Interp *interp,
1700 int argc, CONST84 char *argv[])
1701 {
1702 enum Integrator_type integrator = UNKNOWN;
1703 char buf[MAXIMUM_NUMERIC_LENGTH]; /* string to hold integer */
1704 CONST84 char *engine = NULL;
1705 int result = 0; /* 0 = FALSE; 1 = TRUE */
1706 long i0=(-1), i1=(-1);
1707 int ifound = 0;
1708 int k;
1709 int moststeps=0;
1710 double dt0=0, dtmin=0, dtmax=0;
1711 CONST84 char *cdt0=NULL, *cdtmin=NULL, *cdtmax=NULL, *cmoststeps=NULL,
1712 *ci0=NULL, *ci1=NULL;
1713
1714 (void)cdata; /* stop gcc whine about unused parameter */
1715
1716 k = 1;
1717 while (k < (argc-1)) { /* arguments come in pairs */
1718 if (strcmp(argv[k],"-engine")==0) {
1719 engine = argv[k+1];
1720 k+=2;
1721 continue;
1722 }
1723 if (strcmp(argv[k],"-i1")==0) {
1724 ci1 = argv[k+1];
1725 k+=2;
1726 continue;
1727 }
1728 if (strcmp(argv[k],"-i0")==0) {
1729 ci0 = argv[k+1];
1730 k+=2;
1731 continue;
1732 }
1733 if (strcmp(argv[k],"-moststeps")==0) {
1734 cmoststeps = argv[k+1];
1735 k+=2;
1736 continue;
1737 }
1738 if (strcmp(argv[k],"-dtmax")==0) {
1739 cdtmax = argv[k+1];
1740 k+=2;
1741 continue;
1742 }
1743 if (strcmp(argv[k],"-dtmin")==0) {
1744 cdtmin = argv[k+1];
1745 k+=2;
1746 continue;
1747 }
1748 if (strcmp(argv[k],"-dt0")==0) {
1749 cdt0 = argv[k+1];
1750 k+=2;
1751 continue;
1752 }
1753 Tcl_AppendResult(interp,argv[0],": unrecognized option ",
1754 argv[k],".",SNULL);
1755 return TCL_ERROR;
1756 }
1757
1758 if (engine != NULL && strncmp(engine,"BLSODE",3)==0) {
1759 integrator = BLSODE;
1760 ifound=1;
1761 }
1762 if (!ifound) {
1763 Tcl_SetResult(interp, "Unsupported integrator", TCL_STATIC);
1764 Tcl_AppendResult(interp," ",engine,SNULL);
1765 return TCL_ERROR;
1766 }
1767 if (ci0 != NULL && ci1 != NULL) {
1768 /* get i0, i1 if both supplied. */
1769 long i;
1770 if (Tcl_ExprLong(interp,ci0,&i)==TCL_ERROR|| i<0) {
1771 Tcl_ResetResult(interp);
1772 Tcl_SetResult(interp, "integrate_setup: index i0 invalid", TCL_STATIC);
1773 return TCL_ERROR;
1774 }
1775 i0=i;
1776 if (Tcl_ExprLong(interp,ci1,&i)==TCL_ERROR|| i<i0) {
1777 Tcl_ResetResult(interp);
1778 Tcl_SetResult(interp, "integrate_setup: index i1 invalid", TCL_STATIC);
1779 return TCL_ERROR;
1780 }
1781 i1=i;
1782 }
1783 if (cdt0 != NULL) {
1784 if (Tcl_GetDouble(interp,cdt0,&dt0) != TCL_OK) {
1785 Tcl_ResetResult(interp);
1786 Tcl_AppendResult(interp, "integrate_setup: initial step length invalid",
1787 " (",cdt0,")", SNULL);
1788 return TCL_ERROR;
1789 }
1790 }
1791 if (cdtmin != NULL) {
1792 if (Tcl_GetDouble(interp,cdtmin,&dtmin) != TCL_OK || dtmin < 0) {
1793 Tcl_ResetResult(interp);
1794 Tcl_AppendResult(interp, "integrate_setup: minimum step length invalid",
1795 " (",cdtmin,")", SNULL);
1796 return TCL_ERROR;
1797 }
1798 }
1799 if (cdtmax != NULL) {
1800 if (Tcl_GetDouble(interp,cdtmax,&dtmax) != TCL_OK || dtmax < dtmin) {
1801 Tcl_ResetResult(interp);
1802 Tcl_AppendResult(interp, "integrate_setup: maximum step length invalid",
1803 " (",cdtmax,")", SNULL);
1804 return TCL_ERROR;
1805 }
1806 }
1807 if (cmoststeps != NULL) {
1808 if (Tcl_GetInt(interp,cmoststeps,&moststeps) != TCL_OK || moststeps < 0) {
1809 Tcl_ResetResult(interp);
1810 Tcl_AppendResult(interp, "integrate_setup: maximum internal steps bad",
1811 " (",cmoststeps,")", SNULL);
1812 return TCL_ERROR;
1813 }
1814 }
1815 result = Integ_setup(integrator,i0,i1,dt0,dtmin,dtmax,moststeps);
1816 sprintf(buf, "%d", result);
1817 Tcl_SetResult(interp, buf, TCL_VOLATILE);
1818 return TCL_OK;
1819 }
1820
1821 /*
1822 * Deallocates any memory used and sets all integration global points
1823 * to NULL.
1824 */
1825 static int Integ_Cleanup() {
1826 g_intgsys_cur = NULL;
1827 g_intginst_d_x = NULL;
1828 if(g_intginst_d_y_vars !=NULL) {
1829 ascfree(g_intginst_d_y_vars);
1830 }
1831 if(g_intginst_d_dydx_vars !=NULL) {
1832 ascfree(g_intginst_d_dydx_vars);
1833 }
1834 if(g_intginst_d_obs_vars !=NULL) {
1835 ascfree(g_intginst_d_obs_vars);
1836 }
1837 g_intginst_d_y_vars=NULL;
1838 g_intginst_d_dydx_vars=NULL;
1839 g_intginst_d_obs_vars=NULL;
1840
1841 /*
1842 * Cleanup the derivative stuff.
1843 */
1844 if (g_intg_diff.input_indices) {
1845 ascfree((char *)g_intg_diff.input_indices);
1846 }
1847 if (g_intg_diff.output_indices) {
1848 ascfree((char *)g_intg_diff.output_indices);
1849 }
1850 if (g_intg_diff.y_vars) {
1851 ascfree((char *)g_intg_diff.y_vars);
1852 }
1853 if (g_intg_diff.ydot_vars) {
1854 ascfree((char *)g_intg_diff.ydot_vars);
1855 }
1856 if (g_intg_diff.rlist) {
1857 ascfree((char *)g_intg_diff.rlist);
1858 }
1859 if (g_intg_diff.dydx_dx) {
1860 DestroyDenseMatrix(g_intg_diff.dydx_dx, g_intg_diff.n_eqns);
1861 }
1862
1863 g_intg_diff.input_indices = NULL;
1864 g_intg_diff.output_indices = NULL;
1865 g_intg_diff.y_vars = NULL;
1866 g_intg_diff.ydot_vars = NULL;
1867 g_intg_diff.dydx_dx = NULL;
1868 g_intg_diff.rlist = NULL;
1869 g_intg_diff.n_eqns = 0L;
1870
1871 print_debugstring("integrate_cleanup\n");
1872 return 0;
1873 }
1874
1875 /********************************************************************/
1876
1877 int Asc_IntegCleanupCmd(ClientData cdata,Tcl_Interp *interp,
1878 int argc, CONST84 char *argv[])
1879 {
1880 (void)cdata; /* stop gcc whine about unused parameter */
1881 (void)argv; /* stop gcc whine about unused parameter */
1882
1883 if (argc!=1) {
1884 Tcl_SetResult(interp, "integrate_cleanup takes no arguments", TCL_STATIC);
1885 return TCL_ERROR;
1886 }
1887
1888 Integ_Cleanup();
1889 return TCL_OK;
1890 }

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