/[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 101 - (show annotations) (download) (as text)
Sat Dec 10 04:22:07 2005 UTC (15 years, 4 months ago) by jds
File MIME type: text/x-csrc
File size: 27263 byte(s)
A little more progress killing compiler warnings.
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 = 1.0;
730 UNUSED_PARAMETER(rel);
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 = NULL;
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_assert(ninputs >= 0);
819 Asc_Panic(2, "ExtRel_MapDataToMtx",
820 "ExtRel_MapDataToMtx is totally broken");
821 for (c=0;c<(unsigned long)ninputs;c++) {
822 inst = (struct Instance *)gl_fetch(inputlist,c+1);
823 /*
824 var = var_instance(inst);
825 */
826 used = var_apply_filter(var,d->filter);
827 if (used) {
828 d->nz.col = mtx_org_to_col(d->mtx,var_sindex(var));
829 value = ptr[c] + mtx_value(d->mtx,&(d->nz));
830 mtx_set_value(d->mtx,&(d->nz), value);
831 }
832 }
833 }
834
835
836 /*
837 * ExtRel Finite Differencing.
838 *
839 * This routine actually does the finite differencing.
840 * The jacobian is a single contiguous vector. We load information
841 * in it *row* wise. If we have noutputs x ninputs = 3 x 4, variables,
842 * then jacobian entry 4,
843 * would correspond to jacobian[1][0], i.e., = 0.5 for this eg.
844 *
845 * 2.0 9.0 4.0 6.0
846 * 0.5 1.3 0.0 9.7
847 * 80 7.0 1.0 2.5
848 *
849 * 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
850 * [0][0] [1][0] [2][1]
851 *
852 * When we are finite differencing variable c, we will be loading
853 * jacobian positions c, c+ninputs, c+2*ninputs ....
854 */
855
856 static double CalculateInterval(double varvalue)
857 {
858 (void)varvalue; /* stop gcc whine about unused parameter */
859
860 return (1.0e-05);
861 }
862
863 static int32 ExtRel__FDiff(struct Slv_Interp *slv_interp,
864 int32 (*eval_func) (/* ARGS */),
865 int32 ninputs, int32 noutputs,
866 double *inputs, double *outputs,
867 double *jacobian)
868 {
869 int32 c1,c2, nok = 0;
870 double *tmp_vector;
871 double *ptr;
872 double old_x,interval,value;
873
874 tmp_vector = (double *)asccalloc(noutputs,sizeof(double));
875 for (c1=0;c1<ninputs;c1++) {
876 old_x = inputs[c1]; /* perturb x */
877 interval = CalculateInterval(old_x);
878 inputs[c1] = old_x + interval;
879 nok = (*eval_func)(slv_interp, ninputs, noutputs, /* call routine */
880 inputs, tmp_vector, jacobian);
881 if (nok) break;
882 ptr = &jacobian[c1];
883 for (c2=0;c2<noutputs;c2++) { /* load jacobian */
884 value = (tmp_vector[c2] - outputs[c2])/interval;
885 *ptr = value;
886 ptr += ninputs;
887 }
888 inputs[c1] = old_x;
889 }
890 ascfree((char *)tmp_vector); /* cleanup */
891 return nok;
892 }
893
894
895 static
896 int32 ExtRel_CalcDeriv(struct rel_relation *rel, struct deriv_data *d)
897 {
898 int32 nok = 0;
899 struct Slv_Interp slv_interp;
900 struct ExtRelCache *cache;
901 struct ExternalFunc *efunc;
902 unsigned long whichvar;
903 int32 (*eval_func)();
904 int32 (*deriv_func)();
905
906 assert(rel_extnodeinfo(rel));
907 cache = rel_extcache(rel);
908 whichvar = rel_extwhichvar(rel);
909 efunc = cache->efunc;
910
911 /*
912 * Check and deal with the special case of the first
913 * computation.
914 */
915 if (cache->first_deriv_eval) {
916 cache->newcalc_done = (unsigned)1;
917 cache->first_deriv_eval = (unsigned)0;
918 }
919
920 /*
921 * If a function evaluation was not recently done, then we
922 * can return the results from the cached jacobian.
923 */
924 if (!cache->newcalc_done) {
925 ExtRel_MapDataToMtx(cache->inputlist, whichvar,
926 cache->ninputs, cache->jacobian, d);
927 return 0;
928 }
929
930 /*
931 * If we are here, we have to do the derivative calculation.
932 * The only thing to determine is whether we do analytical
933 * derivatives (deriv_func != NULL) or finite differencing.
934 * In any case init the interpreter.
935 */
936 Init_Slv_Interp(&slv_interp);
937 slv_interp.deriv_eval = (unsigned)1;
938 slv_interp.user_data = cache->user_data;
939 deriv_func = GetDerivFunc(efunc);
940 if (deriv_func) {
941 nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
942 cache->inputs, cache->outputs, cache->jacobian);
943 if (nok) return nok;
944 }
945 else{
946 eval_func = GetValueFunc(efunc);
947 nok = ExtRel__FDiff(&slv_interp, eval_func,
948 cache->ninputs, cache->noutputs,
949 cache->inputs, cache->outputs, cache->jacobian);
950 if (nok) return nok;
951 }
952
953 /*
954 * Cleanup. Ensure that we update the users data, and load
955 * the main matrix with the derivative information.
956 */
957 cache->user_data = slv_interp.user_data; /* save user info */
958 ExtRel_MapDataToMtx(cache->inputlist, whichvar,
959 cache->ninputs, cache->jacobian, d);
960 return 0;
961 }
962
963
964 /*
965 * ExtRel Deriv routines.
966 *
967 * This is the entry point for most cases. ExtRel_CalcDeriv depends
968 * on ExtRel_Evaluate being called immediately before it.
969 */
970
971 real64 ExtRel_Diffs_RHS(struct rel_relation *rel,
972 var_filter_t *filter,
973 int32 row,
974 mtx_matrix_t mtx)
975 {
976 int32 nok;
977 real64 res;
978 struct deriv_data data;
979
980 data.filter = filter;
981 data.mtx = mtx;
982 data.nz.row = row;
983
984 res = ExtRel_Evaluate_RHS(rel);
985 nok = ExtRel_CalcDeriv(rel,&data);
986 if (nok)
987 return 0.0;
988 else
989 return res;
990 }
991
992
993 /*
994 * FIX FIX FIX
995 */
996 real64 ExtRel_Diffs_LHS(struct rel_relation *rel, var_filter_t *filter,
997 int32 row, mtx_matrix_t mtx)
998 {
999 UNUSED_PARAMETER(rel);
1000 UNUSED_PARAMETER(filter);
1001 UNUSED_PARAMETER(row);
1002 UNUSED_PARAMETER(mtx);
1003 return 1.0;
1004 }
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014

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