/[ascend]/trunk/ascend4/solver/rel.c
ViewVC logotype

Annotation of /trunk/ascend4/solver/rel.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 27093 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 /*
2     * SLV: Ascend Numeric Solver
3     * by Karl Michael Westerberg
4     * Created: 2/6/90
5     * Version: $Revision: 1.32 $
6     * Version control file: $RCSfile: rel.c,v $
7     * Date last modified: $Date: 1998/01/29 00:42:28 $
8     * Last modified by: $Author: ballan $
9     *
10     * This file is part of the SLV solver.
11     *
12     * Copyright (C) 1990 Karl Michael Westerberg
13     * Copyright (C) 1993 Joseph Zaher
14     * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15     *
16     * The SLV solver is free software; you can redistribute
17     * it and/or modify it under the terms of the GNU General Public License as
18     * published by the Free Software Foundation; either version 2 of the
19     * License, or (at your option) any later version.
20     *
21     * The SLV solver is distributed in hope that it will be
22     * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24     * General Public License for more details.
25     *
26     * You should have received a copy of the GNU General Public License
27     * along with the program; if not, write to the Free Software Foundation,
28     * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
29     * COPYING. COPYING is found in ../compiler.
30     *
31     */
32    
33     #include <math.h>
34     #include <stdarg.h>
35     #include "utilities/ascConfig.h"
36     #include "utilities/ascPanic.h"
37     #include "utilities/ascMalloc.h"
38     #include "utilities/mem.h"
39     #include "general/list.h"
40     #include "general/dstring.h"
41     #include "compiler/compiler.h"
42     #include "compiler/symtab.h"
43     #include "compiler/fractions.h"
44     #include "compiler/instance_enum.h"
45     #include "compiler/symtab.h"
46     #include "compiler/extfunc.h"
47     #include "compiler/extcall.h"
48     #include "compiler/functype.h"
49     #include "compiler/safe.h"
50     #include "compiler/dimen.h"
51     #include "compiler/types.h"
52     #include "compiler/find.h"
53     #include "compiler/atomvalue.h"
54     #include "compiler/instquery.h"
55     #include "compiler/mathinst.h"
56     #include "compiler/parentchild.h"
57     #include "compiler/instance_io.h"
58     #include "compiler/relation_type.h"
59     #include "compiler/relation.h"
60     #include "compiler/relation_util.h"
61     #include "compiler/packages.h"
62     #define _SLV_SERVER_C_SEEN_ /* for the extrel stuff in header */
63     #include "solver/mtx.h"
64     #include "solver/slv_types.h"
65     #include "solver/var.h"
66     #include "solver/rel.h"
67     #include "solver/discrete.h"
68     #include "solver/conditional.h"
69     #include "solver/logrel.h"
70     #include "solver/bnd.h"
71     #include "solver/slv_server.h"
72    
73     #define IPTR(i) ((struct Instance *)(i))
74     #define REIMPLEMENT 0 /* if set to 1, compiles code tagged with it. */
75     #define REL_DEBUG FALSE
76    
77     /* define symchar names needed */
78     static symchar *g_strings[1];
79     #define INCLUDED_R g_strings[0]
80    
81     static const struct rel_relation g_rel_defaults = {
82     NULL, /* instance */
83     NULL, /* extnode */
84     NULL, /* incidence */
85     e_undefined, /* e_token */
86     0, /* n_incidences */
87     -1, /* mindex */
88     -1, /* sindex */
89     -1, /* model index */
90     (REL_INCLUDED) /* flags */
91     };
92     /*
93     * Don't forget to update the
94     * initialization when the structure
95     * is modified.
96     */
97    
98     static struct rel_relation *copy(const struct rel_relation *rel)
99     {
100     struct rel_relation *newrel;
101     newrel = (struct rel_relation *)ascmalloc( sizeof(struct rel_relation) );
102     *newrel = *rel;
103     return(newrel);
104     }
105    
106     static struct rel_relation *
107     rel_create_extnode(struct rel_relation * rel, struct ExtCallNode *ext)
108     {
109     struct rel_extnode *nodeinfo;
110    
111     nodeinfo = (struct rel_extnode *)ascmalloc(sizeof(struct rel_extnode));
112     nodeinfo->whichvar = (int)ExternalCallVarIndex(ext);
113     nodeinfo->cache = NULL;
114     rel->nodeinfo = nodeinfo;
115     return rel;
116     }
117    
118     struct rel_relation *
119     rel_create(SlvBackendToken instance, struct rel_relation *newrel)
120     {
121     CONST struct relation *instance_relation;
122     struct ExtCallNode *ext;
123     enum Expr_enum ctype;
124    
125     if (newrel==NULL) {
126     newrel = copy(&g_rel_defaults); /* malloc the relation */
127     } else {
128     *newrel = g_rel_defaults; /* init the space we've been sent */
129     }
130     assert(newrel!=NULL);
131     newrel->instance = instance;
132     instance_relation = GetInstanceRelation(IPTR(instance),&ctype);
133     switch (ctype) {
134     case e_token:
135     newrel->type = e_rel_token;
136     break;
137     case e_opcode:
138     Asc_Panic(2, "rel_create", "switching on e_opcode");
139     break;
140     case e_glassbox:
141     newrel->type = e_rel_glassbox;
142     break;
143     case e_blackbox:
144     newrel->type = e_rel_blackbox;
145     ext = BlackBoxExtCall(instance_relation);
146     if (ext) {
147     newrel = rel_create_extnode(newrel,ext);
148     } else {
149     newrel->nodeinfo = NULL;
150     }
151     break;
152     default:
153     FPRINTF(stderr,"Unknown relation type in rel_create\n");
154     break;
155     }
156     return(newrel);
157     }
158    
159     SlvBackendToken rel_instance(struct rel_relation *rel)
160     {
161     if (rel==NULL) return NULL;
162     return (SlvBackendToken) rel->instance;
163     }
164     void rel_write_name(slv_system_t sys,struct rel_relation *rel,FILE *fp)
165     {
166     if (rel == NULL || fp==NULL) return;
167     if (sys!=NULL) {
168     WriteInstanceName(fp,rel_instance(rel),slv_instance(sys));
169     } else {
170     WriteInstanceName(fp,rel_instance(rel),NULL);
171     }
172     }
173    
174     void rel_destroy(struct rel_relation *rel)
175     {
176     struct Instance *inst;
177     switch (rel->type) {
178     case e_rel_token:
179     break;
180     default:
181     break;
182     }
183     if (rel->nodeinfo) {
184     rel->nodeinfo->cache = NULL;
185     ascfree(rel->nodeinfo);
186     rel->nodeinfo = NULL;
187     }
188     ascfree((POINTER)rel->incidence);
189     inst = IPTR(rel->instance);
190     if (inst) {
191     if (GetInterfacePtr(inst)==rel) { /* we own the pointer -- reset it. */
192     SetInterfacePtr(inst,NULL);
193     }
194     }
195     ascfree((POINTER)rel);
196     }
197    
198     uint32 rel_flags( struct rel_relation *rel)
199     {
200     return rel->flags;
201     }
202    
203     void rel_set_flags(struct rel_relation *rel, uint32 flags)
204     {
205     rel->flags = flags;
206     }
207    
208     uint32 rel_flagbit(struct rel_relation *rel, uint32 one)
209     {
210     if (rel==NULL || rel->instance == NULL) {
211     FPRINTF(stderr,"ERROR: rel_flagbit called with bad var.\n");
212     return 0;
213     }
214     return (rel->flags & one);
215     }
216    
217     void rel_set_flagbit(struct rel_relation *rel, uint32 field,uint32 one)
218     {
219     if (one) {
220     rel->flags |= field;
221     } else {
222     rel->flags &= ~field;
223     }
224     }
225    
226     /*
227     * External relations access. See far below for more details.
228     */
229    
230     struct rel_extnode *rel_extnodeinfo( struct rel_relation *rel)
231     {
232     return(rel->nodeinfo);
233     }
234    
235     int32 rel_extwhichvar( struct rel_relation *rel)
236     {
237     if (rel->nodeinfo) {
238     return(rel->nodeinfo->whichvar);
239     } else {
240     return 0; /* never a valid index */
241     }
242     }
243    
244     struct ExtRelCache *rel_extcache( struct rel_relation *rel)
245     {
246     if (rel->nodeinfo) {
247     return(rel->nodeinfo->cache);
248     } else {
249     return NULL;
250     }
251     }
252    
253     void rel_set_extnodeinfo( struct rel_relation *rel, struct rel_extnode *nodeinfo)
254     {
255     rel->nodeinfo = nodeinfo;
256     }
257    
258     void rel_set_extwhichvar(struct rel_relation *rel, int32 whichvar)
259     {
260     rel->nodeinfo->whichvar = whichvar;
261     }
262    
263     void rel_set_extcache( struct rel_relation *rel,struct ExtRelCache * cache)
264     {
265     rel->nodeinfo->cache = cache;
266     }
267    
268     /*
269     * End of External relations access.
270     */
271    
272     /* this needs to be reimplemented properly. */
273     boolean rel_less(struct rel_relation *rel)
274     {
275     switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
276     case e_notequal:
277     case e_less:
278     case e_lesseq:
279     return(TRUE);
280     default:
281     return(FALSE);
282     }
283     }
284    
285     /* this needs to be reimplemented properly. */
286     boolean rel_equal( struct rel_relation *rel)
287     {
288     switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
289     case e_equal:
290     case e_lesseq:
291     case e_greatereq:
292     return(TRUE);
293     default:
294     return(FALSE);
295     }
296     }
297    
298     boolean rel_greater(struct rel_relation *rel)
299     {
300     switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
301     case e_notequal:
302     case e_greater:
303     case e_greatereq:
304     return(TRUE);
305     default:
306     return(FALSE);
307     }
308     }
309    
310     static enum rel_enum compenum2relenum(enum Expr_enum t)
311     {
312     switch (t) {
313     case e_equal:
314     return e_rel_equal;
315     case e_less:
316     return e_rel_less;
317     case e_greater:
318     return e_rel_greater;
319     case e_lesseq:
320     return e_rel_lesseq;
321     case e_greatereq:
322     return e_rel_greatereq;
323     default:
324     FPRINTF(ASCERR,"ERROR (rel.c): compenum2relenum miscalled.\n");
325     return e_rel_not_equal;
326     }
327     }
328     enum rel_enum rel_relop( struct rel_relation *rel)
329     {
330     return
331     compenum2relenum(RelationRelop(
332     GetInstanceRelationOnly(IPTR(rel->instance))));
333     }
334    
335     char *rel_make_name(slv_system_t sys,struct rel_relation *rel)
336     {
337     return WriteInstanceNameString(IPTR(rel->instance),IPTR(slv_instance(sys)));
338     }
339    
340     int32 rel_mindex( struct rel_relation *rel)
341     {
342     return( rel->mindex );
343     }
344    
345     void rel_set_mindex( struct rel_relation *rel, int32 index)
346     {
347     rel->mindex = index;
348     }
349    
350     int32 rel_sindex( const struct rel_relation *rel)
351     {
352     return( rel->sindex );
353     }
354    
355     void rel_set_sindex( struct rel_relation *rel, int32 index)
356     {
357     rel->sindex = index;
358     }
359    
360     int32 rel_model(const struct rel_relation *rel)
361     {
362     return((const int32) rel->model );
363     }
364    
365     void rel_set_model( struct rel_relation *rel, int32 index)
366     {
367     rel->model = index;
368     }
369    
370     real64 rel_residual(struct rel_relation *rel)
371     {
372     return( RelationResidual(GetInstanceRelationOnly(IPTR(rel->instance))));
373     }
374    
375     void rel_set_residual( struct rel_relation *rel, real64 residual)
376     {
377     struct relation *reln;
378     reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
379     SetRelationResidual(reln,residual);
380     }
381    
382     real64 rel_multiplier(struct rel_relation *rel)
383     {
384     return( RelationMultiplier(GetInstanceRelationOnly(IPTR(rel->instance))));
385     }
386    
387     void rel_set_multiplier(struct rel_relation *rel, real64 multiplier)
388     {
389     struct relation *reln;
390     reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
391     SetRelationMultiplier(reln,multiplier);
392     }
393    
394     real64 rel_nominal( struct rel_relation *rel)
395     {
396     return( RelationNominal(GetInstanceRelationOnly(IPTR(rel->instance))));
397     }
398    
399     void rel_set_nominal( struct rel_relation *rel, real64 nominal)
400     {
401     struct relation *reln;
402     reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
403     SetRelationNominal(reln,nominal);
404     }
405    
406     /* to bad there's no entry point that rel must call before being used
407     * generally, like the FindType checking stuff in var.c
408     */
409     static void check_included_flag(void){
410     if (INCLUDED_R == NULL || AscFindSymbol(INCLUDED_R) == NULL) {
411     INCLUDED_R = AddSymbolL("included",8);
412     }
413     }
414     uint32 rel_included( struct rel_relation *rel)
415     {
416     struct Instance *c;
417     check_included_flag();
418     c = ChildByChar(IPTR(rel->instance),INCLUDED_R);
419     if( c == NULL ) {
420     FPRINTF(stderr,"ERROR: (rel) rel_included\n");
421     FPRINTF(stderr," No 'included' field in relation.\n");
422     WriteInstance(stderr,IPTR(rel->instance));
423     return FALSE;
424     }
425     rel_set_flagbit(rel,REL_INCLUDED,GetBooleanAtomValue(c));
426     return( GetBooleanAtomValue(c) );
427     }
428    
429     void rel_set_included( struct rel_relation *rel, uint32 included)
430     {
431     struct Instance *c;
432     check_included_flag();
433     c = ChildByChar(IPTR(rel->instance),INCLUDED_R);
434     if( c == NULL ) {
435     FPRINTF(stderr,"ERROR: (rel) rel_set_included\n");
436     FPRINTF(stderr," No 'included' field in relation.\n");
437     WriteInstance(stderr,IPTR(rel->instance));
438     return;
439     }
440     SetBooleanAtomValue(c,included,(unsigned)0);
441     rel_set_flagbit(rel,REL_INCLUDED,included);
442     }
443    
444     int32 rel_apply_filter( const struct rel_relation *rel, rel_filter_t *filter)
445     {
446     if (rel==NULL || filter==NULL) {
447     FPRINTF(stderr,"rel_apply_filter miscalled with NULL\n");
448     return FALSE;
449     }
450     /* AND to mask off irrelevant bits in flags and match value, then compare */
451     return ( (filter->matchbits & rel->flags) ==
452     (filter->matchbits & filter->matchvalue)
453     );
454     }
455    
456     int32 rel_n_incidencesF(struct rel_relation *rel)
457     {
458     if (rel!=NULL) return rel->n_incidences;
459     FPRINTF(stderr,"rel_n_incidences miscalled with NULL\n");
460     return 0;
461     }
462     void rel_set_incidencesF(struct rel_relation *rel,int32 n,
463     struct var_variable **i)
464     {
465     if(rel!=NULL && n >=0) {
466     if (n && i==NULL) {
467     FPRINTF(stderr,"rel_set_incidence miscalled with NULL ilist\n");
468     }
469     rel->n_incidences = n;
470     rel->incidence = i;
471     return;
472     }
473     FPRINTF(stderr,"rel_set_incidence miscalled with NULL or n < 0\n");
474     }
475     const struct var_variable **rel_incidence_list( struct rel_relation *rel)
476     {
477     if (rel==NULL) return NULL;
478     return( (const struct var_variable **)rel->incidence );
479     }
480    
481     struct var_variable **rel_incidence_list_to_modify( struct rel_relation *rel)
482     {
483     if (rel==NULL) return NULL;
484     return( (struct var_variable **)rel->incidence );
485     }
486    
487     #if KILL
488     expr_t rel_lhs( struct rel_relation *rel)
489     {
490     if (rel==NULL) return NULL;
491     return( rel->lhs );
492     }
493    
494     expr_t rel_rhs( struct rel_relation *rel)
495     {
496     if (rel==NULL) return NULL;
497     return( rel->rhs );
498     }
499     #endif /* KILL */
500    
501     /*
502     * External Relations Cache for solvers.
503     * by Kirk Andre Abbott
504     * Created: 08/10/94
505     */
506    
507     double g_external_tolerance = 1.0e-12;
508    
509     struct ExtRelCache *CreateExtRelCache(struct ExtCallNode *ext)
510     {
511     struct ExtRelCache *result;
512     unsigned long n_input_args, n_output_args;
513     int32 ninputs, noutputs;
514    
515     assert(ext!=NULL);
516     result = (struct ExtRelCache *)ascmalloc(sizeof(struct ExtRelCache));
517     result->nodestamp = ExternalCallNodeStamp(ext);
518     result->efunc = ExternalCallExtFunc(ext);
519     result->data = ExternalCallDataInstance(ext);
520     result->arglist = ExternalCallArgList(ext);
521     result->user_data = NULL; /* absolutely important to be
522     * initialized to NULL ! */
523    
524     n_input_args = NumberInputArgs(result->efunc);
525     n_output_args = NumberOutputArgs(result->efunc);
526    
527     /*
528     * We own the result of the LinearizeArgList call.
529     */
530     result->inputlist = LinearizeArgList(result->arglist,1,n_input_args);
531     ninputs = (int)gl_length(result->inputlist);
532     noutputs = (int)CountNumberOfArgs(result->arglist,n_input_args+1,
533     n_input_args+n_output_args);
534     result->ninputs = ninputs;
535     result->noutputs = noutputs;
536    
537     result->inputs = (double *)asccalloc(ninputs,sizeof(double));
538     result->outputs = (double *)asccalloc(noutputs,sizeof(double));
539     result->jacobian = (double *)asccalloc(ninputs*noutputs,sizeof(double));
540     /*
541     * Setup default flags for controlling calculations.
542     */
543     result->newcalc_done = (unsigned)1;
544     return result;
545     }
546    
547    
548     struct ExtRelCache *CreateCacheFromInstance(SlvBackendToken relinst)
549     {
550     struct ExtCallNode *ext;
551     struct ExtRelCache *cache;
552     CONST struct relation *reln;
553     enum Expr_enum type;
554    
555     assert(relinst != NULL && InstanceKind(IPTR(relinst))==REL_INST);
556     reln = GetInstanceRelation(IPTR(relinst),&type);
557     if (type!=e_blackbox) {
558     FPRINTF(stderr,"Invalid relation type in (CreateCacheFromInstance)\n");
559     return NULL;
560     }
561     ext = BlackBoxExtCall(reln);
562     cache = CreateExtRelCache(ext);
563     return cache;
564     }
565    
566     void ExtRel_DestroyCache(struct ExtRelCache *cache)
567     {
568     if (cache) {
569     cache->nodestamp=-1;
570     cache->efunc = NULL;
571     cache->data = NULL;
572     gl_destroy(cache->inputlist); /* we own the list */
573     ascfree(cache->inputs);
574     ascfree(cache->outputs);
575     ascfree(cache->jacobian);
576     cache->inputlist = NULL;
577     cache->inputs = NULL;
578     cache->outputs = NULL;
579     cache->jacobian = NULL;
580     ascfree(cache);
581     }
582     }
583    
584    
585     /*
586     * ExtRel_PreSolve:
587     *
588     * To deal with the first time we also want to do arguement
589     * checking, and then turn off the first_func_eval flag.
590     * Turn on the newcalc_done flag. The rationale behind this is
591     * as follows:
592     * The solver at the moment does not treat an external relation
593     * specially, i.e., as a block. It also calls for its functions
594     * a relation at a time. However the external relations compute
595     * all their outputs at once. So as not to do unnecessary
596     * recalculations we use these flag bits. We set newcalc_done
597     * initially to true, so as to force *at least* one calculation
598     * of the external relations. By similar reasoning first_func_eval (done)
599     * is set to false.
600     */
601     int32 ExtRel_PreSolve(struct ExtRelCache *cache, int32 setup)
602     {
603     struct ExternalFunc *efunc;
604     struct Slv_Interp slv_interp;
605    
606     int32 (*init_func)(struct Slv_Interp *,
607     struct Instance *,struct gl_list_t *);
608     int32 nok = 0;
609    
610     if (!cache) return 1;
611     efunc = cache->efunc;
612     init_func = GetInitFunc(efunc);
613     Init_Slv_Interp(&slv_interp);
614     slv_interp.nodestamp = cache->nodestamp;
615     slv_interp.user_data = cache->user_data;
616     if (setup) {
617     slv_interp.first_call = (unsigned)1;
618     slv_interp.last_call = (unsigned)0;
619     slv_interp.check_args = (unsigned)1;
620     }
621     else{
622     slv_interp.first_call = (unsigned)0;
623     slv_interp.last_call = (unsigned)1;
624     slv_interp.check_args = (unsigned)0;
625     }
626     nok = (*init_func)(&slv_interp,cache->data,cache->arglist);
627     if (nok)
628     return 1;
629    
630     /*
631     * Save the user's data and update our status.
632     */
633     cache->user_data = slv_interp.user_data;
634     cache->newcalc_done = (unsigned)1; /* force at least one calculation */
635     cache->first_func_eval = (unsigned)0;
636     return 0;
637     }
638    
639    
640     static int32 ArgsDifferent(double new, double old)
641     {
642     if (fabs(new - old) >= g_external_tolerance) {
643     return 1;
644     } else {
645     return 0;
646     }
647     }
648    
649     real64 ExtRel_Evaluate_RHS(struct rel_relation *rel)
650     {
651     struct Slv_Interp slv_interp;
652     struct ExtRelCache *cache;
653     struct ExternalFunc *efunc;
654     struct Instance *arg;
655     struct gl_list_t *inputlist;
656     double value;
657     int32 c,ninputs;
658     int32 nok;
659     unsigned long whichvar;
660     int32 newcalc_reqd=0;
661    
662     /* badly need to use a function prototype typedef in the header. */
663     int32 (*eval_func)(struct Slv_Interp *,
664     int, int,
665     double *, double *, double *);
666    
667     assert(rel_extnodeinfo(rel));
668     cache = rel_extcache(rel);
669     efunc = cache->efunc;
670     eval_func = GetValueFunc(efunc);
671     inputlist = cache->inputlist;
672     ninputs = cache->ninputs;
673     whichvar = (unsigned long)rel_extwhichvar(rel);
674    
675     /*
676     * The determination of whether a new calculation is required needs
677     * some more thought. One thing we should insist upon is that all
678     * the relations for an external relation are forced into the same
679     * block.
680     */
681     for (c=0;c<ninputs;c++) {
682     arg = (struct Instance *)gl_fetch(inputlist,c+1);
683     value = RealAtomValue(arg);
684     if (ArgsDifferent(value,cache->inputs[c])) {
685     newcalc_reqd = 1;
686     cache->inputs[c] = value;
687     }
688     }
689     value = 0;
690     /*
691     * Do the calculations if necessary. Note that we have to *ensure*
692     * that we send the user the information that he provided to us.
693     * We have to update our user_data after each call of the user's function
694     * as he might move information around (not smart but possible), on us.
695     * If a function call is made, mark a new calculation as having been,
696     * done, otherwise dont.
697     */
698     if (newcalc_reqd) {
699     Init_Slv_Interp(&slv_interp);
700     slv_interp.nodestamp = cache->nodestamp;
701     slv_interp.user_data = cache->user_data;
702     slv_interp.func_eval = (unsigned)1;
703    
704     nok = (*eval_func)(&slv_interp, ninputs, cache->noutputs,
705     cache->inputs, cache->outputs, cache->jacobian);
706     value = cache->outputs[whichvar - ninputs - 1];
707     cache->newcalc_done = (unsigned)1; /* newcalc done */
708     cache->user_data = slv_interp.user_data; /* update user_data */
709     }
710     else{
711     value = cache->outputs[whichvar - ninputs - 1];
712     cache->newcalc_done = (unsigned)0; /* a result was simply returned */
713     }
714    
715     #ifdef KAA_DEBUG
716     FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_RHS result ->%g\n",
717     result);
718     #endif /* KAA_DEBUG */
719    
720     return value;
721     }
722    
723    
724     /*
725     * FIX FIX FIX
726     */
727     real64 ExtRel_Evaluate_LHS(struct rel_relation *rel)
728     {
729     real64 res;
730     res = 1.0;
731     FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_LHS result ->%g\n",
732     res);
733     return res;
734     }
735    
736    
737     /*
738     * ExtRel_Gradient evaluation routines.
739     *
740     * The following code implements gradient evaluation routines for
741     * external relations. The routines here will prepare the arguements
742     * and call a user supplied derivative routine, is same is non-NULL.
743     * If it is the user supplied function evaluation routine will be
744     * used to compute the gradients via finite differencing.
745     * The current solver will not necessarily call for the derivative
746     * all at once. This makes it necessary to do the gradient
747     * computations (user supplied or finite difference), and to cache
748     * away the results. Based on calculation flags, the appropriate
749     * *row* of this cached jacobian will be extracted and mapped to the
750     * main solve matrix.
751     *
752     * The cached jacobian is a contiguous vector ninputs*noutputs long
753     * and is loaded row wise. Indexing starts from 0. Each row corresponds
754     * to the partial derivatives of the output variable (associated with
755     * that row, wrt to all its incident input variables.
756     *
757     * Careful attention needs to be paid to the way this jacobian is
758     * loaded/unloaded, because of the multiple indexing schemes in use.
759     * i.e, arglist's and inputlists index 1..nvars, whereas all numeric
760     * vectors number from 0.
761     */
762    
763     struct deriv_data {
764     var_filter_t *filter;
765     mtx_matrix_t mtx;
766     mtx_coord_t nz;
767     };
768    
769    
770     /*
771     * ExtRel_MapDataToMtx
772     *
773     * This function attempts to map the information from the contiguous
774     * jacobian back into the main matrix, based on the current row <=>
775     * whichvar. The jacobian assumed to have been calculated.
776     * Since we are operating a relation at a time, we have to find
777     * out where to index into our jacobian. This index is computed as
778     * follows:
779     *
780     * index = (whichvar - ninputs - 1) * ninputs
781     *
782     * Example: a problem with 4 inputs, 3 outputs and whichvar = 6.
783     * with the counting for vars 1..nvars, but the jacobian indexing
784     * starting from 0 (c-wise).
785     *
786     * V-------- first output variable
787     * 1 2 3 4 5 6 7
788     * ^---------------- whichvar
789     * ------------------- grads for whichvar = 6
790     * | | | |
791     * v v v v
792     * index = 0 1 2 3 4 5 6 7 8 9 10 11
793     * jacobian = 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
794     *
795     * Hence jacobian index = (6 - 4 - 1) * 4 = 4
796     *
797     *
798     * THIS FUNCTION IS TOTALLY AND COMPLETELY BROKEN.
799     */
800    
801     static void ExtRel_MapDataToMtx(struct gl_list_t *inputlist,
802     unsigned long whichvar,
803     int32 ninputs,
804     double *jacobian,
805     struct deriv_data *d)
806     {
807     struct Instance *inst;
808     struct var_variable *var;
809     double value, *ptr;
810     boolean used;
811     unsigned long c;
812     int32 index;
813    
814     index = ((int)whichvar - ninputs - 1) * ninputs;
815     ptr = &jacobian[index];
816    
817     /* this is totally broken, thanks to kirk making the var=instance assumption */
818     Asc_Panic(2, "ExtRel_MapDataToMtx",
819     "ExtRel_MapDataToMtx is totally broken");
820     for (c=0;c<ninputs;c++) {
821     inst = (struct Instance *)gl_fetch(inputlist,c+1);
822     /*
823     var = var_instance(inst);
824     */
825     used = var_apply_filter(var,d->filter);
826     if (used) {
827     d->nz.col = mtx_org_to_col(d->mtx,var_sindex(var));
828     value = ptr[c] + mtx_value(d->mtx,&(d->nz));
829     mtx_set_value(d->mtx,&(d->nz), value);
830     }
831     }
832     }
833    
834    
835     /*
836     * ExtRel Finite Differencing.
837     *
838     * This routine actually does the finite differencing.
839     * The jacobian is a single contiguous vector. We load information
840     * in it *row* wise. If we have noutputs x ninputs = 3 x 4, variables,
841     * then jacobian entry 4,
842     * would correspond to jacobian[1][0], i.e., = 0.5 for this eg.
843     *
844     * 2.0 9.0 4.0 6.0
845     * 0.5 1.3 0.0 9.7
846     * 80 7.0 1.0 2.5
847     *
848     * 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
849     * [0][0] [1][0] [2][1]
850     *
851     * When we are finite differencing variable c, we will be loading
852     * jacobian positions c, c+ninputs, c+2*ninputs ....
853     */
854    
855     static double CalculateInterval(double varvalue)
856     {
857     (void)varvalue; /* stop gcc whine about unused parameter */
858    
859     return (1.0e-05);
860     }
861    
862     static int32 ExtRel__FDiff(struct Slv_Interp *slv_interp,
863     int32 (*eval_func) (/* ARGS */),
864     int32 ninputs, int32 noutputs,
865     double *inputs, double *outputs,
866     double *jacobian)
867     {
868     int32 c1,c2, nok = 0;
869     double *tmp_vector;
870     double *ptr;
871     double old_x,interval,value;
872    
873     tmp_vector = (double *)asccalloc(noutputs,sizeof(double));
874     for (c1=0;c1<ninputs;c1++) {
875     old_x = inputs[c1]; /* perturb x */
876     interval = CalculateInterval(old_x);
877     inputs[c1] = old_x + interval;
878     nok = (*eval_func)(slv_interp, ninputs, noutputs, /* call routine */
879     inputs, tmp_vector, jacobian);
880     if (nok) break;
881     ptr = &jacobian[c1];
882     for (c2=0;c2<noutputs;c2++) { /* load jacobian */
883     value = (tmp_vector[c2] - outputs[c2])/interval;
884     *ptr = value;
885     ptr += ninputs;
886     }
887     inputs[c1] = old_x;
888     }
889     ascfree((char *)tmp_vector); /* cleanup */
890     return nok;
891     }
892    
893    
894     static
895     int32 ExtRel_CalcDeriv(struct rel_relation *rel, struct deriv_data *d)
896     {
897     int32 nok = 0;
898     struct Slv_Interp slv_interp;
899     struct ExtRelCache *cache;
900     struct ExternalFunc *efunc;
901     unsigned long whichvar;
902     int32 (*eval_func)();
903     int32 (*deriv_func)();
904    
905     assert(rel_extnodeinfo(rel));
906     cache = rel_extcache(rel);
907     whichvar = rel_extwhichvar(rel);
908     efunc = cache->efunc;
909    
910     /*
911     * Check and deal with the special case of the first
912     * computation.
913     */
914     if (cache->first_deriv_eval) {
915     cache->newcalc_done = (unsigned)1;
916     cache->first_deriv_eval = (unsigned)0;
917     }
918    
919     /*
920     * If a function evaluation was not recently done, then we
921     * can return the results from the cached jacobian.
922     */
923     if (!cache->newcalc_done) {
924     ExtRel_MapDataToMtx(cache->inputlist, whichvar,
925     cache->ninputs, cache->jacobian, d);
926     return 0;
927     }
928    
929     /*
930     * If we are here, we have to do the derivative calculation.
931     * The only thing to determine is whether we do analytical
932     * derivatives (deriv_func != NULL) or finite differencing.
933     * In any case init the interpreter.
934     */
935     Init_Slv_Interp(&slv_interp);
936     slv_interp.deriv_eval = (unsigned)1;
937     slv_interp.user_data = cache->user_data;
938     deriv_func = GetDerivFunc(efunc);
939     if (deriv_func) {
940     nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
941     cache->inputs, cache->outputs, cache->jacobian);
942     if (nok) return nok;
943     }
944     else{
945     eval_func = GetValueFunc(efunc);
946     nok = ExtRel__FDiff(&slv_interp, eval_func,
947     cache->ninputs, cache->noutputs,
948     cache->inputs, cache->outputs, cache->jacobian);
949     if (nok) return nok;
950     }
951    
952     /*
953     * Cleanup. Ensure that we update the users data, and load
954     * the main matrix with the derivative information.
955     */
956     cache->user_data = slv_interp.user_data; /* save user info */
957     ExtRel_MapDataToMtx(cache->inputlist, whichvar,
958     cache->ninputs, cache->jacobian, d);
959     return 0;
960     }
961    
962    
963     /*
964     * ExtRel Deriv routines.
965     *
966     * This is the entry point for most cases. ExtRel_CalcDeriv depends
967     * on ExtRel_Evaluate being called immediately before it.
968     */
969    
970     real64 ExtRel_Diffs_RHS(struct rel_relation *rel,
971     var_filter_t *filter,
972     int32 row,
973     mtx_matrix_t mtx)
974     {
975     int32 nok;
976     real64 res;
977     struct deriv_data data;
978    
979     data.filter = filter;
980     data.mtx = mtx;
981     data.nz.row = row;
982    
983     res = ExtRel_Evaluate_RHS(rel);
984     nok = ExtRel_CalcDeriv(rel,&data);
985     if (nok)
986     return 0.0;
987     else
988     return res;
989     }
990    
991    
992     /*
993     * FIX FIX FIX
994     */
995     real64 ExtRel_Diffs_LHS(struct rel_relation *rel, var_filter_t *filter,
996     int32 row, mtx_matrix_t mtx)
997     {
998     return 1.0;
999     }
1000    
1001    
1002    
1003    
1004    
1005    
1006    
1007    
1008    
1009    

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