/[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 694 - (show annotations) (download) (as text)
Thu Jun 22 07:54:53 2006 UTC (17 years, 11 months ago) by johnpye
File MIME type: text/x-csrc
File size: 27512 byte(s)
Lots of tidying up to relation_util code. (Most commments but a 
little refactoring of repeated code chunks.)
Added 'IMPORT' to the Gedit syntax definition.
Added error_reporters in place of some console output.
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/expr_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 "mtx.h"
64 #include "slv_types.h"
65 #include "var.h"
66 #include "rel.h"
67 #include "discrete.h"
68 #include "conditional.h"
69 #include "logrel.h"
70 #include "bnd.h"
71 #include "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 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"blackbox support not yet implemented");
145 newrel->type = e_rel_blackbox;
146 ext = BlackBoxExtCall(instance_relation);
147 if (ext) {
148 newrel = rel_create_extnode(newrel,ext);
149 } else {
150 newrel->nodeinfo = NULL;
151 }
152 break;
153 default:
154 FPRINTF(stderr,"Unknown relation type in rel_create\n");
155 break;
156 }
157 return(newrel);
158 }
159
160 SlvBackendToken rel_instance(struct rel_relation *rel)
161 {
162 if (rel==NULL) return NULL;
163 return (SlvBackendToken) rel->instance;
164 }
165 void rel_write_name(slv_system_t sys,struct rel_relation *rel,FILE *fp)
166 {
167 if (rel == NULL || fp==NULL) return;
168 if (sys!=NULL) {
169 WriteInstanceName(fp,rel_instance(rel),slv_instance(sys));
170 } else {
171 WriteInstanceName(fp,rel_instance(rel),NULL);
172 }
173 }
174
175 void rel_destroy(struct rel_relation *rel)
176 {
177 struct Instance *inst;
178 switch (rel->type) {
179 case e_rel_token:
180 break;
181 default:
182 break;
183 }
184 if (rel->nodeinfo) {
185 rel->nodeinfo->cache = NULL;
186 ascfree(rel->nodeinfo);
187 rel->nodeinfo = NULL;
188 }
189 ascfree((POINTER)rel->incidence);
190 inst = IPTR(rel->instance);
191 if (inst) {
192 if (GetInterfacePtr(inst)==rel) { /* we own the pointer -- reset it. */
193 SetInterfacePtr(inst,NULL);
194 }
195 }
196 ascfree((POINTER)rel);
197 }
198
199 uint32 rel_flags( struct rel_relation *rel)
200 {
201 return rel->flags;
202 }
203
204 void rel_set_flags(struct rel_relation *rel, uint32 flags)
205 {
206 rel->flags = flags;
207 }
208
209 uint32 rel_flagbit(struct rel_relation *rel, uint32 one)
210 {
211 if (rel==NULL || rel->instance == NULL) {
212 FPRINTF(stderr,"ERROR: rel_flagbit called with bad var.\n");
213 return 0;
214 }
215 return (rel->flags & one);
216 }
217
218 void rel_set_flagbit(struct rel_relation *rel, uint32 field,uint32 one)
219 {
220 if (one) {
221 rel->flags |= field;
222 } else {
223 rel->flags &= ~field;
224 }
225 }
226
227 /*
228 * External relations access. See far below for more details.
229 */
230
231 struct rel_extnode *rel_extnodeinfo( struct rel_relation *rel)
232 {
233 return(rel->nodeinfo);
234 }
235
236 int32 rel_extwhichvar( struct rel_relation *rel)
237 {
238 if (rel->nodeinfo) {
239 return(rel->nodeinfo->whichvar);
240 } else {
241 return 0; /* never a valid index */
242 }
243 }
244
245 struct ExtRelCache *rel_extcache( struct rel_relation *rel)
246 {
247 if (rel->nodeinfo) {
248 return(rel->nodeinfo->cache);
249 } else {
250 return NULL;
251 }
252 }
253
254 void rel_set_extnodeinfo( struct rel_relation *rel, struct rel_extnode *nodeinfo)
255 {
256 rel->nodeinfo = nodeinfo;
257 }
258
259 void rel_set_extwhichvar(struct rel_relation *rel, int32 whichvar)
260 {
261 rel->nodeinfo->whichvar = whichvar;
262 }
263
264 void rel_set_extcache( struct rel_relation *rel,struct ExtRelCache * cache)
265 {
266 rel->nodeinfo->cache = cache;
267 }
268
269 /*
270 * End of External relations access.
271 */
272
273 /* this needs to be reimplemented properly. */
274 boolean rel_less(struct rel_relation *rel)
275 {
276 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
277 case e_notequal:
278 case e_less:
279 case e_lesseq:
280 return(TRUE);
281 default:
282 return(FALSE);
283 }
284 }
285
286 /* this needs to be reimplemented properly. */
287 boolean rel_equal( struct rel_relation *rel)
288 {
289 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
290 case e_equal:
291 case e_lesseq:
292 case e_greatereq:
293 return(TRUE);
294 default:
295 return(FALSE);
296 }
297 }
298
299 boolean rel_greater(struct rel_relation *rel)
300 {
301 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
302 case e_notequal:
303 case e_greater:
304 case e_greatereq:
305 return(TRUE);
306 default:
307 return(FALSE);
308 }
309 }
310
311 static enum rel_enum compenum2relenum(enum Expr_enum t)
312 {
313 switch (t) {
314 case e_equal:
315 return e_rel_equal;
316 case e_less:
317 return e_rel_less;
318 case e_greater:
319 return e_rel_greater;
320 case e_lesseq:
321 return e_rel_lesseq;
322 case e_greatereq:
323 return e_rel_greatereq;
324 default:
325 FPRINTF(ASCERR,"ERROR (rel.c): compenum2relenum miscalled.\n");
326 return e_rel_not_equal;
327 }
328 }
329 enum rel_enum rel_relop( struct rel_relation *rel)
330 {
331 return
332 compenum2relenum(RelationRelop(
333 GetInstanceRelationOnly(IPTR(rel->instance))));
334 }
335
336 char *rel_make_name(slv_system_t sys,struct rel_relation *rel)
337 {
338 return WriteInstanceNameString(IPTR(rel->instance),IPTR(slv_instance(sys)));
339 }
340
341 int32 rel_mindex( struct rel_relation *rel)
342 {
343 return( rel->mindex );
344 }
345
346 void rel_set_mindex( struct rel_relation *rel, int32 index)
347 {
348 rel->mindex = index;
349 }
350
351 int32 rel_sindex( const struct rel_relation *rel)
352 {
353 return( rel->sindex );
354 }
355
356 void rel_set_sindex( struct rel_relation *rel, int32 index)
357 {
358 rel->sindex = index;
359 }
360
361 int32 rel_model(const struct rel_relation *rel)
362 {
363 return((const int32) rel->model );
364 }
365
366 void rel_set_model( struct rel_relation *rel, int32 index)
367 {
368 rel->model = index;
369 }
370
371 real64 rel_residual(struct rel_relation *rel)
372 {
373 return( RelationResidual(GetInstanceRelationOnly(IPTR(rel->instance))));
374 }
375
376 void rel_set_residual( struct rel_relation *rel, real64 residual)
377 {
378 struct relation *reln;
379 reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
380 SetRelationResidual(reln,residual);
381 }
382
383 real64 rel_multiplier(struct rel_relation *rel)
384 {
385 return( RelationMultiplier(GetInstanceRelationOnly(IPTR(rel->instance))));
386 }
387
388 void rel_set_multiplier(struct rel_relation *rel, real64 multiplier)
389 {
390 struct relation *reln;
391 reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
392 SetRelationMultiplier(reln,multiplier);
393 }
394
395 real64 rel_nominal( struct rel_relation *rel)
396 {
397 return( RelationNominal(GetInstanceRelationOnly(IPTR(rel->instance))));
398 }
399
400 void rel_set_nominal( struct rel_relation *rel, real64 nominal)
401 {
402 struct relation *reln;
403 reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
404 SetRelationNominal(reln,nominal);
405 }
406
407 /* to bad there's no entry point that rel must call before being used
408 * generally, like the FindType checking stuff in var.c
409 */
410 static void check_included_flag(void){
411 if (INCLUDED_R == NULL || AscFindSymbol(INCLUDED_R) == NULL) {
412 INCLUDED_R = AddSymbolL("included",8);
413 }
414 }
415 uint32 rel_included( struct rel_relation *rel)
416 {
417 struct Instance *c;
418 check_included_flag();
419 c = ChildByChar(IPTR(rel->instance),INCLUDED_R);
420 if( c == NULL ) {
421 FPRINTF(stderr,"ERROR: (rel) rel_included\n");
422 FPRINTF(stderr," No 'included' field in relation.\n");
423 WriteInstance(stderr,IPTR(rel->instance));
424 return FALSE;
425 }
426 rel_set_flagbit(rel,REL_INCLUDED,GetBooleanAtomValue(c));
427 return( GetBooleanAtomValue(c) );
428 }
429
430 void rel_set_included( struct rel_relation *rel, uint32 included)
431 {
432 struct Instance *c;
433 check_included_flag();
434 c = ChildByChar(IPTR(rel->instance),INCLUDED_R);
435 if( c == NULL ) {
436 FPRINTF(stderr,"ERROR: (rel) rel_set_included\n");
437 FPRINTF(stderr," No 'included' field in relation.\n");
438 WriteInstance(stderr,IPTR(rel->instance));
439 return;
440 }
441 SetBooleanAtomValue(c,included,(unsigned)0);
442 rel_set_flagbit(rel,REL_INCLUDED,included);
443 }
444
445 int32 rel_apply_filter( const struct rel_relation *rel, rel_filter_t *filter)
446 {
447 if (rel==NULL || filter==NULL) {
448 FPRINTF(stderr,"rel_apply_filter miscalled with NULL\n");
449 return FALSE;
450 }
451 /* AND to mask off irrelevant bits in flags and match value, then compare */
452 return ( (filter->matchbits & rel->flags) ==
453 (filter->matchbits & filter->matchvalue)
454 );
455 }
456
457 int32 rel_n_incidencesF(struct rel_relation *rel)
458 {
459 if (rel!=NULL) return rel->n_incidences;
460 FPRINTF(stderr,"rel_n_incidences miscalled with NULL\n");
461 return 0;
462 }
463 void rel_set_incidencesF(struct rel_relation *rel,int32 n,
464 struct var_variable **i)
465 {
466 if(rel!=NULL && n >=0) {
467 if (n && i==NULL) {
468 FPRINTF(stderr,"rel_set_incidence miscalled with NULL ilist\n");
469 }
470 rel->n_incidences = n;
471 rel->incidence = i;
472 return;
473 }
474 FPRINTF(stderr,"rel_set_incidence miscalled with NULL or n < 0\n");
475 }
476 const struct var_variable **rel_incidence_list( struct rel_relation *rel)
477 {
478 if (rel==NULL) return NULL;
479 return( (const struct var_variable **)rel->incidence );
480 }
481
482 struct var_variable **rel_incidence_list_to_modify( struct rel_relation *rel)
483 {
484 if (rel==NULL) return NULL;
485 return( (struct var_variable **)rel->incidence );
486 }
487
488 #if KILL
489 expr_t rel_lhs( struct rel_relation *rel)
490 {
491 if (rel==NULL) return NULL;
492 return( rel->lhs );
493 }
494
495 expr_t rel_rhs( struct rel_relation *rel)
496 {
497 if (rel==NULL) return NULL;
498 return( rel->rhs );
499 }
500 #endif /* KILL */
501
502 /*
503 * External Relations Cache for solvers.
504 * by Kirk Andre Abbott
505 * Created: 08/10/94
506 */
507
508 double g_external_tolerance = 1.0e-12;
509
510 struct ExtRelCache *CreateExtRelCache(struct ExtCallNode *ext)
511 {
512 struct ExtRelCache *result;
513 unsigned long n_input_args, n_output_args;
514 int32 ninputs, noutputs;
515
516 assert(ext!=NULL);
517 result = (struct ExtRelCache *)ascmalloc(sizeof(struct ExtRelCache));
518 result->nodestamp = ExternalCallNodeStamp(ext);
519 result->efunc = ExternalCallExtFunc(ext);
520 result->data = ExternalCallDataInstance(ext);
521 result->arglist = ExternalCallArgList(ext);
522 result->user_data = NULL; /* absolutely important to be
523 * initialized to NULL ! */
524
525 n_input_args = NumberInputArgs(result->efunc);
526 n_output_args = NumberOutputArgs(result->efunc);
527
528 /*
529 * We own the result of the LinearizeArgList call.
530 */
531 result->inputlist = LinearizeArgList(result->arglist,1,n_input_args);
532 ninputs = (int)gl_length(result->inputlist);
533 noutputs = (int)CountNumberOfArgs(result->arglist,n_input_args+1,
534 n_input_args+n_output_args);
535 result->ninputs = ninputs;
536 result->noutputs = noutputs;
537
538 result->inputs = ASC_NEW_ARRAY_CLEAR(double,ninputs);
539 result->outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);
540 result->jacobian = (double *)asccalloc(ninputs*noutputs,sizeof(double));
541 /*
542 * Setup default flags for controlling calculations.
543 */
544 result->newcalc_done = (unsigned)1;
545 return result;
546 }
547
548
549 struct ExtRelCache *CreateCacheFromInstance(SlvBackendToken relinst)
550 {
551 struct ExtCallNode *ext;
552 struct ExtRelCache *cache;
553 CONST struct relation *reln;
554 enum Expr_enum type;
555
556 assert(relinst != NULL && InstanceKind(IPTR(relinst))==REL_INST);
557 reln = GetInstanceRelation(IPTR(relinst),&type);
558 if (type!=e_blackbox) {
559 FPRINTF(stderr,"Invalid relation type in (CreateCacheFromInstance)\n");
560 return NULL;
561 }
562 ext = BlackBoxExtCall(reln);
563 cache = CreateExtRelCache(ext);
564 return cache;
565 }
566
567 void ExtRel_DestroyCache(struct ExtRelCache *cache)
568 {
569 if (cache) {
570 cache->nodestamp=-1;
571 cache->efunc = NULL;
572 cache->data = NULL;
573 gl_destroy(cache->inputlist); /* we own the list */
574 ascfree(cache->inputs);
575 ascfree(cache->outputs);
576 ascfree(cache->jacobian);
577 cache->inputlist = NULL;
578 cache->inputs = NULL;
579 cache->outputs = NULL;
580 cache->jacobian = NULL;
581 ascfree(cache);
582 }
583 }
584
585
586 /*
587 * ExtRel_PreSolve:
588 *
589 * To deal with the first time we also want to do arguement
590 * checking, and then turn off the first_func_eval flag.
591 * Turn on the newcalc_done flag. The rationale behind this is
592 * as follows:
593 * The solver at the moment does not treat an external relation
594 * specially, i.e., as a block. It also calls for its functions
595 * a relation at a time. However the external relations compute
596 * all their outputs at once. So as not to do unnecessary
597 * recalculations we use these flag bits. We set newcalc_done
598 * initially to true, so as to force *at least* one calculation
599 * of the external relations. By similar reasoning first_func_eval (done)
600 * is set to false.
601 */
602 int32 ExtRel_PreSolve(struct ExtRelCache *cache, int32 setup)
603 {
604 struct ExternalFunc *efunc;
605 struct Slv_Interp slv_interp;
606
607 ExtBBoxInitFunc *init_func;
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 /*
618 slv_interp.first_call = (unsigned)1;
619 slv_interp.last_call = (unsigned)0;
620 slv_interp.check_args = (unsigned)1;
621 */
622 slv_interp.task = bb_first_call;
623 }
624 else{
625 slv_interp.task = bb_last_call;
626 /*
627 slv_interp.first_call = (unsigned)0;
628 slv_interp.last_call = (unsigned)1;
629 slv_interp.check_args = (unsigned)0;
630 */
631 }
632 nok = (*init_func)(&slv_interp,cache->data,cache->arglist);
633 if (nok) {
634 return 1;
635 }
636
637 /*
638 * Save the user's data and update our status.
639 */
640 cache->user_data = slv_interp.user_data;
641 cache->newcalc_done = (unsigned)1; /* force at least one calculation */
642 cache->first_func_eval = (unsigned)0;
643 return 0;
644 }
645
646
647 static int32 ArgsDifferent(double new, double old)
648 {
649 if (fabs(new - old) >= g_external_tolerance) {
650 return 1;
651 } else {
652 return 0;
653 }
654 }
655
656 real64 ExtRel_Evaluate_RHS(struct rel_relation *rel)
657 {
658 struct Slv_Interp slv_interp;
659 struct ExtRelCache *cache;
660 struct ExternalFunc *efunc;
661 struct Instance *arg;
662 struct gl_list_t *inputlist;
663 double value;
664 int32 c,ninputs;
665 int32 nok;
666 unsigned long whichvar;
667 int32 newcalc_reqd=0;
668
669 /* badly need to use a function prototype typedef in the header. */
670 int32 (*eval_func)(struct Slv_Interp *,
671 int, int,
672 double *, double *, double *);
673
674 assert(rel_extnodeinfo(rel));
675 cache = rel_extcache(rel);
676 efunc = cache->efunc;
677 eval_func = GetValueFunc(efunc);
678 inputlist = cache->inputlist;
679 ninputs = cache->ninputs;
680 whichvar = (unsigned long)rel_extwhichvar(rel);
681
682 /*
683 * The determination of whether a new calculation is required needs
684 * some more thought. One thing we should insist upon is that all
685 * the relations for an external relation are forced into the same
686 * block.
687 */
688 for (c=0;c<ninputs;c++) {
689 arg = (struct Instance *)gl_fetch(inputlist,c+1);
690 value = RealAtomValue(arg);
691 if (ArgsDifferent(value,cache->inputs[c])) {
692 newcalc_reqd = 1;
693 cache->inputs[c] = value;
694 }
695 }
696 value = 0;
697 /*
698 * Do the calculations if necessary. Note that we have to *ensure*
699 * that we send the user the information that he provided to us.
700 * We have to update our user_data after each call of the user's function
701 * as he might move information around (not smart but possible), on us.
702 * If a function call is made, mark a new calculation as having been,
703 * done, otherwise dont.
704 */
705 if (newcalc_reqd) {
706 Init_Slv_Interp(&slv_interp);
707 slv_interp.nodestamp = cache->nodestamp;
708 slv_interp.user_data = cache->user_data;
709 /*
710 slv_interp.func_eval = (unsigned)1;
711 */
712 slv_interp.task = bb_func_eval;
713
714 nok = (*eval_func)(&slv_interp, ninputs, cache->noutputs,
715 cache->inputs, cache->outputs, cache->jacobian);
716 value = cache->outputs[whichvar - ninputs - 1];
717 cache->newcalc_done = (unsigned)1; /* newcalc done */
718 cache->user_data = slv_interp.user_data; /* update user_data */
719 }
720 else{
721 value = cache->outputs[whichvar - ninputs - 1];
722 cache->newcalc_done = (unsigned)0; /* a result was simply returned */
723 }
724
725 #ifdef KAA_DEBUG
726 FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_RHS result ->%g\n",
727 result);
728 #endif /* KAA_DEBUG */
729
730 return value;
731 }
732
733
734 /*
735 * FIX FIX FIX
736 */
737 real64 ExtRel_Evaluate_LHS(struct rel_relation *rel)
738 {
739 real64 res = 0.0;
740 /*
741 conceptual equation is 0 = f(y[i], x[]);
742 where f is the residual function for the ith output y[i] and all input x[j]
743 */
744 UNUSED_PARAMETER(rel);
745 FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_LHS result ->%g\n",
746 res);
747 return res;
748 }
749
750
751 /*
752 * ExtRel_Gradient evaluation routines.
753 *
754 * The following code implements gradient evaluation routines for
755 * external relations. The routines here will prepare the arguements
756 * and call a user supplied derivative routine, is same is non-NULL.
757 * If it is the user supplied function evaluation routine will be
758 * used to compute the gradients via finite differencing.
759 * The current solver will not necessarily call for the derivative
760 * all at once. This makes it necessary to do the gradient
761 * computations (user supplied or finite difference), and to cache
762 * away the results. Based on calculation flags, the appropriate
763 * *row* of this cached jacobian will be extracted and mapped to the
764 * main solve matrix.
765 *
766 * The cached jacobian is a contiguous vector ninputs*noutputs long
767 * and is loaded row wise. Indexing starts from 0. Each row corresponds
768 * to the partial derivatives of the output variable (associated with
769 * that row, wrt to all its incident input variables.
770 *
771 * Careful attention needs to be paid to the way this jacobian is
772 * loaded/unloaded, because of the multiple indexing schemes in use.
773 * i.e, arglist's and inputlists index 1..nvars, whereas all numeric
774 * vectors number from 0.
775 */
776
777 struct deriv_data {
778 var_filter_t *filter;
779 mtx_matrix_t mtx;
780 mtx_coord_t nz;
781 };
782
783
784 /*
785 * ExtRel_MapDataToMtx
786 *
787 * This function attempts to map the information from the contiguous
788 * jacobian back into the main matrix, based on the current row <=>
789 * whichvar. The jacobian assumed to have been calculated.
790 * Since we are operating a relation at a time, we have to find
791 * out where to index into our jacobian. This index is computed as
792 * follows:
793 *
794 * index = (whichvar - ninputs - 1) * ninputs
795 *
796 * Example: a problem with 4 inputs, 3 outputs and whichvar = 6.
797 * with the counting for vars 1..nvars, but the jacobian indexing
798 * starting from 0 (c-wise).
799 *
800 * V-------- first output variable
801 * 1 2 3 4 5 6 7
802 * ^---------------- whichvar
803 * ------------------- grads for whichvar = 6
804 * | | | |
805 * v v v v
806 * index = 0 1 2 3 4 5 6 7 8 9 10 11
807 * jacobian = 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
808 *
809 * Hence jacobian index = (6 - 4 - 1) * 4 = 4
810 *
811 *
812 * THIS FUNCTION IS TOTALLY AND COMPLETELY BROKEN.
813 */
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 {
872 (void)varvalue; /* stop gcc whine about unused parameter */
873
874 return (1.0e-05);
875 }
876
877 static int32 ExtRel__FDiff(struct Slv_Interp *slv_interp,
878 int32 (*eval_func) (/* ARGS */),
879 int32 ninputs, int32 noutputs,
880 double *inputs, double *outputs,
881 double *jacobian)
882 {
883 int32 c1,c2, nok = 0;
884 double *tmp_vector;
885 double *ptr;
886 double old_x,interval,value;
887
888 tmp_vector = ASC_NEW_ARRAY_CLEAR(double,noutputs);
889 for (c1=0;c1<ninputs;c1++) {
890 old_x = inputs[c1]; /* perturb x */
891 interval = CalculateInterval(old_x);
892 inputs[c1] = old_x + interval;
893 nok = (*eval_func)(slv_interp, ninputs, noutputs, /* call routine */
894 inputs, tmp_vector, jacobian);
895 if (nok) break;
896 ptr = &jacobian[c1];
897 for (c2=0;c2<noutputs;c2++) { /* load jacobian */
898 value = (tmp_vector[c2] - outputs[c2])/interval;
899 *ptr = value;
900 ptr += ninputs;
901 }
902 inputs[c1] = old_x;
903 }
904 ascfree((char *)tmp_vector); /* cleanup */
905 return nok;
906 }
907
908
909 static
910 int32 ExtRel_CalcDeriv(struct rel_relation *rel, struct deriv_data *d)
911 {
912 int32 nok = 0;
913 struct Slv_Interp slv_interp;
914 struct ExtRelCache *cache;
915 struct ExternalFunc *efunc;
916 unsigned long whichvar;
917 int32 (*eval_func)();
918 int32 (*deriv_func)();
919
920 assert(rel_extnodeinfo(rel));
921 cache = rel_extcache(rel);
922 whichvar = rel_extwhichvar(rel);
923 efunc = cache->efunc;
924
925 /*
926 * Check and deal with the special case of the first
927 * computation.
928 */
929 if (cache->first_deriv_eval) {
930 cache->newcalc_done = (unsigned)1;
931 cache->first_deriv_eval = (unsigned)0;
932 }
933
934 /*
935 * If a function evaluation was not recently done, then we
936 * can return the results from the cached jacobian.
937 */
938 if (!cache->newcalc_done) {
939 ExtRel_MapDataToMtx(cache->inputlist, whichvar,
940 cache->ninputs, cache->jacobian, d);
941 return 0;
942 }
943
944 /*
945 * If we are here, we have to do the derivative calculation.
946 * The only thing to determine is whether we do analytical
947 * derivatives (deriv_func != NULL) or finite differencing.
948 * In any case init the interpreter.
949 */
950 Init_Slv_Interp(&slv_interp);
951 /*
952 slv_interp.deriv_eval = (unsigned)1;
953 */
954 slv_interp.task = bb_deriv_eval;
955 slv_interp.user_data = cache->user_data;
956 deriv_func = GetDerivFunc(efunc);
957 if (deriv_func) {
958 nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
959 cache->inputs, cache->outputs, cache->jacobian);
960 if (nok) return nok;
961 }
962 else{
963 eval_func = GetValueFunc(efunc);
964 nok = ExtRel__FDiff(&slv_interp, eval_func,
965 cache->ninputs, cache->noutputs,
966 cache->inputs, cache->outputs, cache->jacobian);
967 if (nok) return nok;
968 }
969
970 /*
971 * Cleanup. Ensure that we update the users data, and load
972 * the main matrix with the derivative information.
973 */
974 cache->user_data = slv_interp.user_data; /* save user info */
975 ExtRel_MapDataToMtx(cache->inputlist, whichvar,
976 cache->ninputs, cache->jacobian, d);
977 return 0;
978 }
979
980
981 /*
982 * ExtRel Deriv routines.
983 *
984 * This is the entry point for most cases. ExtRel_CalcDeriv depends
985 * on ExtRel_Evaluate being called immediately before it.
986 */
987
988 real64 ExtRel_Diffs_RHS(struct rel_relation *rel,
989 var_filter_t *filter,
990 int32 row,
991 mtx_matrix_t mtx)
992 {
993 int32 nok;
994 real64 res;
995 struct deriv_data data;
996
997 data.filter = filter;
998 data.mtx = mtx;
999 data.nz.row = row;
1000
1001 res = ExtRel_Evaluate_RHS(rel);
1002 nok = ExtRel_CalcDeriv(rel,&data);
1003 if (nok)
1004 return 0.0;
1005 else
1006 return res;
1007 }
1008
1009
1010 /*
1011 * FIX FIX FIX
1012 */
1013 real64 ExtRel_Diffs_LHS(struct rel_relation *rel, var_filter_t *filter,
1014 int32 row, mtx_matrix_t mtx)
1015 {
1016 UNUSED_PARAMETER(rel);
1017 UNUSED_PARAMETER(filter);
1018 UNUSED_PARAMETER(row);
1019 UNUSED_PARAMETER(mtx);
1020 return 1.0;
1021 }
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031

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