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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 571 - (hide 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 johnpye 571 /*
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