/[ascend]/trunk/base/generic/solver/rel.c
ViewVC logotype

Contents of /trunk/base/generic/solver/rel.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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