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

Contents of /trunk/base/generic/solver/relman.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: 26484 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 /*
2 * SLV: Ascend Nonlinear Solver
3 * by Karl Michael Westerberg
4 * Created: 2/6/90
5 * Version: $Revision: 1.44 $
6 * Version control file: $RCSfile: relman.c,v $
7 * Date last modified: $Date: 1998/04/23 23:56:22 $
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 along with
27 * the program; if not, write to the Free Software Foundation, Inc., 675
28 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
29 * COPYING is found in ../compiler.
30 */
31
32 #include <math.h>
33 #include <utilities/ascConfig.h>
34 #include <compiler/instance_enum.h>
35 #include <compiler/fractions.h>
36 #include <compiler/compiler.h>
37 #include <utilities/ascMalloc.h>
38 #include <compiler/functype.h>
39 #include <compiler/safe.h>
40 #include <general/list.h>
41 #include <compiler/extfunc.h>
42 #include <compiler/dimen.h>
43 #include <compiler/expr_types.h>
44 #include <compiler/find.h>
45 #include <compiler/atomvalue.h>
46 #include <compiler/mathinst.h>
47 #include <compiler/relation_type.h>
48 #include <compiler/relation.h>
49 #include <compiler/relation_util.h>
50 #include <compiler/relation_io.h>
51 #define _SLV_SERVER_C_SEEN_
52 #include "mtx.h"
53 #include "slv_types.h"
54 #include "var.h"
55 #include "rel.h"
56 #include "discrete.h"
57 #include "conditional.h"
58 #include "logrel.h"
59 #include "bnd.h"
60 #include "relman.h"
61 #include <compiler/rootfind.h>
62 #include "slv_server.h"
63 #include <general/mathmacros.h>
64
65 #define IPTR(i) ((struct Instance *)(i))
66
67 #define KILL 0 /* compile dead code if kill = 1 */
68 #define REIMPLEMENT 0 /* code that needs to be reimplemented */
69
70 static POINTER rel_tmpalloc( int nbytes)
71 /**
72 *** Temporarily allocates a given number of bytes. The memory need
73 *** not be freed, but the next call to this function will reuse the
74 *** previous allocation. Memory returned will NOT be zeroed.
75 *** Calling with nbytes==0 will free any memory allocated.
76 **/
77 {
78 static char *ptr = NULL;
79 static int cap = 0;
80
81 if (nbytes) {
82 if( nbytes > cap ) {
83 if( ptr != NULL ) ascfree(ptr);
84 ptr = ASC_NEW_ARRAY(char,nbytes);
85 cap = nbytes;
86 }
87 } else {
88 if (ptr) ascfree(ptr);
89 ptr=NULL;
90 cap=0;
91 }
92 if ( cap >0) {
93 return(ptr);
94 } else {
95 return NULL;
96 }
97 }
98
99 #define rel_tmpalloc_array(nelts,type) \
100 ((nelts) > 0 ? (type *)tmpalloc((nelts)*sizeof(type)) : NULL)
101 /**
102 *** Creates an array of "nelts" objects, each with type "type".
103 **/
104
105 void relman_free_reused_mem(void)
106 {
107 /* rel_tmpalloc(0); */
108 RelationFindRoots(NULL,0,0,0,0,NULL,NULL,NULL);
109 }
110
111
112 #if REIMPLEMENT
113 boolean relman_is_linear( struct rel_relation *rel, var_filter_t *filter)
114 {
115 return (
116 exprman_is_linear(rel,rel_lhs(rel),filter) &&
117 exprman_is_linear(rel,rel_rhs(rel),filter)
118 );
119 }
120
121 real64 relman_linear_coef(struct rel_relation *rel, struct var_variable *var,
122 var_filter_t *filter)
123 {
124 return(
125 exprman_linear_coef(rel,rel_lhs(rel),var,filter) -
126 exprman_linear_coef(rel,rel_rhs(rel),var,filter)
127 );
128 }
129 #endif
130
131 #if KILL
132 void relman_decide_incidence( struct rel_relation *rel)
133 {
134 struct var_variable **list;
135 int c;
136
137 list = rel_incidence_list(rel);
138 for( c=0; list[c] != NULL; c++ )
139 var_set_incident(list[c],TRUE);
140 }
141 #endif
142
143 void relman_get_incidence(struct rel_relation *rel, var_filter_t *filter,
144 mtx_matrix_t mtx)
145 {
146 const struct var_variable **list;
147 mtx_coord_t nz;
148 int c,len;
149
150 assert(rel!=NULL && filter !=NULL && mtx != NULL);
151 nz.row = rel_sindex(rel);
152 len = rel_n_incidences(rel);
153
154 list = rel_incidence_list(rel);
155 for (c=0; c < len; c++) {
156 if( var_apply_filter(list[c],filter) ) {
157 nz.col = var_sindex(list[c]);
158 mtx_fill_org_value(mtx,&nz,1.0);
159 }
160 }
161 }
162
163 /*
164 *********************************************************************
165 * Code to deal with glassbox relations processing.
166 *********************************************************************
167 */
168
169 static double dsolve_scratch = 0.0; /* some workspace */
170 #define DSOLVE_TOLERANCE 1.0e-08 /* no longer needed */
171
172 static
173 real64 *relman_glassbox_dsolve(struct rel_relation *rel,
174 struct var_variable *solvefor,
175 int *able,
176 int *nsolns,
177 real64 tolerance)
178 {
179 int n,m,j,dummy = 0;
180 int mode, result;
181 int index,solve_for_index = -1;
182 struct Instance *var;
183 CONST struct relation *cmplr_reln;
184 enum Expr_enum type;
185 CONST struct gl_list_t *vlist;
186 double *f, *x, *g, value;
187 double lower, upper, nominal;
188 ExtEvalFunc **table;
189
190 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
191 assert(type==e_glassbox);
192
193 vlist = RelationVarList(cmplr_reln);
194 m = rel_sindex(rel);
195 n = (int)gl_length(vlist);
196 f = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n));
197 x = &f[1];
198 g = &f[n+1];
199
200 /*
201 * Load x vector, and get the index into the x[vector]
202 * of the variable that we are solving for.
203 */
204 for (j=0;j<n;j++) {
205 var = (struct Instance *)gl_fetch(vlist,(unsigned long)(j+1));
206 x[j] = RealAtomValue(var);
207 if (var_instance(solvefor)== var)
208 solve_for_index = j;
209 }
210
211 /*
212 * Set up bounds and tolerance.
213 */
214 lower = var_lower_bound(solvefor);
215 upper = var_upper_bound(solvefor);
216 nominal = var_nominal(solvefor);
217 /*tolerance = DSOLVE_TOLERANCE; now passed in. */
218
219 /*
220 * Get the evaluation function.
221 */
222 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
223 index = GlassBoxRelIndex(cmplr_reln);
224
225 /** @TODO FIXME
226 Call the rootfinder. A note of CAUTION:
227 zbrent is going to call the evaluation function.
228 It expects, that the results of this evaluation will be
229 written to f[m]. GlassBox relations however always
230 write to slot f[0]. To keep the world happy we send in
231 a dummy = 0. This needs to be made cleaner.
232 */
233 value = zbrent(table[index], &lower, &upper,
234 &mode, &dummy, &solve_for_index,
235 x, x , f, g, &tolerance, &result);
236 if (result==0) {
237 dsolve_scratch = value;
238 ascfree((char *)f);
239 *able = TRUE;
240 *nsolns = 1;
241 return &dsolve_scratch;
242 }
243 else{
244 dsolve_scratch = 0.0;
245 ascfree((char *)f);
246 *able = FALSE;
247 *nsolns = 0;
248 return NULL;
249 }
250 }
251
252
253 #ifdef THIS_IS_AN_UNUSED_FUNCTION
254 static
255 real64 relman_glassbox_eval(struct rel_relation *rel)
256 {
257 int n,m,mode,result;
258 int index,j;
259 CONST struct relation *cmplr_reln;
260 enum Expr_enum type;
261 CONST struct gl_list_t *vlist;
262 double *f, *x, *g, value;
263 ExtEvalFunc **table, *eval;
264
265 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
266 assert(type==e_glassbox);
267
268 vlist = RelationVarList(cmplr_reln);
269 m = rel_sindex(rel);
270 n = (int)gl_length(vlist);
271 /* this needs to be reused memory, fergossake */
272 f = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n)); /* resid */
273 x = &f[1]; /* var values */
274 g = &f[n+1]; /* gradient */
275
276 for (j=0;j<n;j++) {
277 x[j] = RealAtomValue(gl_fetch(vlist,(unsigned long)(j+1)));
278 }
279
280 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
281 index = GlassBoxRelIndex(cmplr_reln);
282
283 /*
284 * Remember to set up the fpe traps etc ....
285 * and to monitor the return flag.
286 */
287 mode = 0;
288 eval = table[index];
289 result = (*eval)(&mode,&m,&n,x,NULL,f,g);
290 value = f[0];
291 ascfree((char *)f);
292 return value;
293 }
294 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
295
296
297 #define BROKENKIRK 0
298 #if BROKENKIRK
299 /* fills filter passing gradient elements to matrix */
300 /* this needs to be buried on the compiler side. */
301 void relman_map_grad2mtx( struct rel_relation *rel, var_filter_t *filter,
302 mtx_matrix_t mtx, CONST struct gl_list_t *varlist, double *g, int m, int n)
303 {
304 mtx_coord_t nz;
305 struct var_variable *var;
306 double value;
307 int j;
308
309 nz.row = m; /* org row of glassbox reln? rel_sindex? */
310
311 for (j=0;j<n;j++) {
312 /* var = (struct Instance *)gl_fetch(varlist,(unsigned long)(j+1)); */
313 var = rel_incidence(rel)[j];
314 if (var_apply_filter(var)) {
315 nz.col = var_sindex(var);
316 value = g[j] /* + mtx_value(mtx,&nz) no longer incrementing row */;
317 mtx_fill_org_value(mtx,&nz, value);
318 }
319 }
320 }
321
322 real64 relman_glassbox_diffs( struct rel_relation *rel,
323 var_filter_t *filter, mtx_matrix_t mtx)
324 {
325 int n,m,mode,result;
326 int index,j;
327 struct Instance *var;
328 CONST struct relation *cmplr_reln;
329 enum Expr_enum type;
330 CONST struct gl_list_t *vlist;
331 double *f, *x, *g, value;
332 ExtEvalFunc **table, *eval;
333
334 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
335 vlist = RelationVarList(cmplr_reln);
336 m = rel_sindex(rel);
337 n = (int)gl_length(vlist);
338 /* this needs to be reused memory ! */
339 f = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n));
340 x = &f[1];
341 g = &f[n+1];
342
343 for (j=0;j<n;j++) {
344 x[j] = RealAtomValue(gl_fetch(vlist,(unsigned long)j+1));
345 }
346
347 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
348 index = GlassBoxRelIndex(cmplr_reln);
349 eval = table[index];
350 result = (*eval)(0,&m,&n,x,NULL,f,g);
351
352 table = GetDerivJumpTable(GlassBoxExtFunc(cmplr_reln));
353
354 mode = 0;
355 eval = table[index];
356 result += (*eval)(&mode,&m,&n,x,NULL,f,g);
357 relman_map_grad2mtx(rel,filter,mtx,vlist,g,m,n);
358
359 /*
360 * Remember to set up the fpe traps etc ....
361 * and to monitor the return flag.
362 */
363 value = f[0];
364 ascfree((char *)f);
365 return value;
366 }
367
368 #else
369 #define relman_glassbox_diffs(rel,filter,mtx) abort()
370 #endif
371
372 /* returns residual; sets *calc_ok = 1 on success */
373 real64 relman_eval(struct rel_relation *rel, int32 *calc_ok, int safe){
374 real64 res;
375 assert(calc_ok!=NULL && rel!=NULL);
376
377 /*
378 token relations
379 */
380 if( rel->type == e_rel_token ){
381 if(!RelationCalcResidualBinary(
382 GetInstanceRelationOnly(IPTR(rel->instance))
383 ,&res)
384 ){
385 *calc_ok = 1;
386 rel_set_residual(rel,res);
387 return res;
388 }
389 }
390
391 if(rel->nodeinfo){
392 CONSOLE_DEBUG("ABOUT TO CALL EXTREL_EVALUATE_RESIDUAL");
393 CONSOLE_DEBUG("REL_RELATION = %p",rel);
394 *calc_ok = 1;
395 res = ExtRel_Evaluate_Residual(rel);
396 return res;
397 }
398
399 /*
400 other types of relations (which include...). This latter approach is
401 apparently older and "reasonably correct, if slow".
402 */
403 CONSOLE_DEBUG("EVALUATE REL_RELATION = %p",rel);
404 if(safe){
405 *calc_ok = (int32)RelationCalcResidualSafe(rel_instance(rel),&res);
406 safe_error_to_stderr( (enum safe_err *)calc_ok );
407
408 /* always set the relation residual when using safe functions */
409 rel_set_residual(rel,res);
410 }else{
411 *calc_ok = RelationCalcResidual(rel_instance(rel),&res);
412
413 /* for unsafe functions, set a fallback value in case of error */
414 if(*calc_ok){
415 res = 1.0e8;
416 }else{
417 rel_set_residual(rel,res);
418 }
419 }
420
421 /* flip the status flag: all values other than safe_ok become 0 */
422 *calc_ok = !(*calc_ok);
423 return res;
424 }
425
426 int32 relman_obj_direction(struct rel_relation *rel)
427 {
428 assert(rel!=NULL);
429 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
430 case e_minimize:
431 return -1;
432 case e_maximize:
433 return 1;
434 default:
435 return 0;
436 }
437 }
438
439 real64 relman_scale(struct rel_relation *rel)
440 {
441 real64 relnom;
442 assert(rel!=NULL);
443 relnom = CalcRelationNominal(rel_instance(rel));
444 if(relnom < 0.0000001) {
445 /* take care of small relnoms and relnom = 0 error returns */
446 relnom = 1.0;
447 }
448 rel_set_nominal(rel,relnom);
449 return relnom;
450 }
451
452 #if REIMPLEMENT /* compiler */
453 real64 relman_diff(struct rel_relation *rel, struct var_variable *var,
454 int safe)
455 {
456 /* FIX FIX FIX meaning kirk couldn't be botghered... */
457 real64 res = 0.0;
458 switch(rel->type) {
459 case e_glassbox:
460 case e_blackbox:
461 FPRINTF(stderr,"relman_diff not yet supported for black and glassbox\n");
462 return 1.0e08;
463 }
464 res -= exprman_diff(rel,rel_rhs(rel),var);
465 res += exprman_diff(rel,rel_lhs(rel),var);
466 return( res );
467 }
468 #endif
469
470 int relman_diff2(struct rel_relation *rel, var_filter_t *filter,
471 real64 *derivatives, int32 *variables,
472 int32 *count, int32 safe)
473 {
474 const struct var_variable **vlist=NULL;
475 real64 *gradient;
476 int32 len,c;
477 int status;
478
479 assert(rel!=NULL && filter!=NULL);
480 len = rel_n_incidences(rel);
481 vlist = rel_incidence_list(rel);
482
483 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
484 assert(gradient !=NULL);
485 *count = 0;
486 if( safe ) {
487 status =(int32)RelationCalcGradientSafe(rel_instance(rel),gradient);
488 safe_error_to_stderr( (enum safe_err *)&status );
489 /* always map when using safe functions */
490 for (c=0; c < len; c++) {
491 if (var_apply_filter(vlist[c],filter)) {
492 variables[*count] = var_sindex(vlist[c]);
493 derivatives[*count] = gradient[c];
494 (*count)++;
495 }
496 }
497 }
498 else {
499 if((status=RelationCalcGradient(rel_instance(rel),gradient)) == 0) {
500 /* successful */
501 for (c=0; c < len; c++) {
502 if (var_apply_filter(vlist[c],filter)) {
503 variables[*count] = var_sindex(vlist[c]);
504 derivatives[*count] = gradient[c];
505 (*count)++;
506 }
507 }
508 }
509 }
510
511 /* flip the status flag */
512 return !status;
513 }
514
515 int relman_diff_grad(struct rel_relation *rel, var_filter_t *filter,
516 real64 *derivatives, int32 *variables_master,
517 int32 *variables_solver, int32 *count, real64 *resid,
518 int32 safe)
519 {
520 const struct var_variable **vlist=NULL;
521 real64 *gradient;
522 int32 len,c;
523 int status;
524
525 assert(rel!=NULL && filter!=NULL);
526 len = rel_n_incidences(rel);
527 vlist = rel_incidence_list(rel);
528
529 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
530 assert(gradient !=NULL);
531 *count = 0;
532 if( safe ) {
533 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),
534 resid,gradient);
535 safe_error_to_stderr( (enum safe_err *)&status );
536 /* always map when using safe functions */
537 for (c=0; c < len; c++) {
538 if (var_apply_filter(vlist[c],filter)) {
539 variables_master[*count] = var_mindex(vlist[c]);
540 variables_solver[*count] = var_sindex(vlist[c]);
541 derivatives[*count] = gradient[c];
542 (*count)++;
543 }
544 }
545 }
546 else {
547 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient))== 0) {
548 /* successful */
549 for (c=0; c < len; c++) {
550 if (var_apply_filter(vlist[c],filter)) {
551 variables_master[*count] = var_mindex(vlist[c]);
552 variables_solver[*count] = var_sindex(vlist[c]);
553 derivatives[*count] = gradient[c];
554 (*count)++;
555 }
556 }
557 }
558 }
559
560 return !status; /* flip the status flag */
561 }
562
563 int32 relman_diff_harwell(struct rel_relation **rlist,
564 var_filter_t *vfilter, rel_filter_t *rfilter,
565 int32 rlen, int32 bias, int32 mORs,
566 real64 *avec, int32 *ivec, int32 *jvec)
567 {
568 const struct var_variable **vlist = NULL;
569 struct rel_relation *rel;
570 real64 residual, *resid;
571 real64 *gradient;
572 mtx_coord_t coord;
573 int32 len,c,r,k;
574 int32 errcnt;
575 enum safe_err status;
576
577 if (rlist == NULL || vfilter == NULL || rfilter == NULL || avec == NULL ||
578 rlen < 0 || mORs >3 || mORs < 0 || bias <0 || bias > 1) {
579 return 1;
580 }
581 resid = &residual;
582 errcnt = k = 0;
583
584 if ( ivec == NULL || jvec == NULL ) {
585 /*_skip stuffing ivec,jvec */
586 for (r=0; r < rlen; r++) {
587 rel = rlist[r];
588 len = rel_n_incidences(rel);
589 vlist = rel_incidence_list(rel);
590 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
591 if (gradient == NULL) {
592 return 1;
593 }
594 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
595 safe_error_to_stderr(&status);
596 if (status) {
597 errcnt--;
598 }
599 for (c=0; c < len; c++) {
600 if (var_apply_filter(vlist[c],vfilter)) {
601 avec[k] = gradient[c];
602 k++;
603 }
604 }
605 }
606 } else {
607 for (r=0; r < rlen; r++) {
608 rel = rlist[r];
609 len = rel_n_incidences(rel);
610 vlist = rel_incidence_list(rel);
611 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
612 if (gradient == NULL) {
613 return 1;
614 }
615 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
616 safe_error_to_stderr(&status);
617 if (status) {
618 errcnt--;
619 }
620 if (mORs & 2) {
621 coord.row = rel_sindex(rel);
622 } else {
623 coord.row = rel_mindex(rel);
624 }
625 for (c=0; c < len; c++) {
626 if (var_apply_filter(vlist[c],vfilter)) {
627 if (mORs & 1) {
628 coord.col = var_sindex(vlist[c]);
629 } else {
630 coord.col = var_mindex(vlist[c]);
631 }
632 avec[k] = gradient[c];
633 ivec[k] = coord.row;
634 jvec[k] = coord.col;
635 k++;
636 }
637 }
638 }
639 }
640 return errcnt;
641 }
642
643 int32 relman_jacobian_count(struct rel_relation **rlist, int32 rlen,
644 var_filter_t *vfilter,
645 rel_filter_t *rfilter, int32 *max)
646 {
647 int32 len, result=0, row, count, c;
648 const struct var_variable **list;
649 struct rel_relation *rel;
650 *max = 0;
651
652 for( row = 0; row < rlen; row++ ) {
653 rel = rlist[row];
654 if (rel_apply_filter(rel,rfilter)) {
655 len = rel_n_incidences(rel);
656 list = rel_incidence_list(rel);
657 count = 0;
658 for (c=0; c < len; c++) {
659 if( var_apply_filter(list[c],vfilter) ) {
660 count++;
661 }
662 }
663 result += count;
664 *max = MAX(*max,count);
665 }
666 }
667 return result;
668 }
669
670 int relman_diffs(struct rel_relation *rel, var_filter_t *filter,
671 mtx_matrix_t mtx, real64 *resid, int safe)
672 {
673 const struct var_variable **vlist=NULL;
674 real64 *gradient;
675 int32 len,c;
676 mtx_coord_t coord;
677 int status;
678
679 assert(rel!=NULL && filter!=NULL && mtx != NULL);
680 len = rel_n_incidences(rel);
681 vlist = rel_incidence_list(rel);
682 coord.row = rel_sindex(rel);
683 assert(coord.row>=0 && coord.row < mtx_order(mtx));
684
685 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
686 assert(gradient !=NULL);
687 if( safe ) {
688 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
689 safe_error_to_stderr( (enum safe_err *)&status );
690 /* always map when using safe functions */
691 for (c=0; c < len; c++) {
692 if (var_apply_filter(vlist[c],filter)) {
693 coord.col = var_sindex(vlist[c]);
694 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
695 mtx_fill_org_value(mtx,&coord,gradient[c]);
696 }
697 }
698 }
699 else {
700 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient)) == 0) {
701 /* successful */
702 for (c=0; c < len; c++) {
703 if (var_apply_filter(vlist[c],filter)) {
704 coord.col = var_sindex(vlist[c]);
705 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
706 mtx_fill_org_value(mtx,&coord,gradient[c]);
707 }
708 }
709 }
710 }
711
712 /* flip the status flag */
713 return !status;
714
715 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
716 switch (rel->type) {
717 case e_token:
718 vlist = rel_incidence_list(rel);
719 res -= exprman_diffs_alt(rel,rel_rhs(rel),filter,row,mtx,vlist);
720 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
721 res += exprman_diffs_alt(rel,rel_lhs(rel),filter,row,mtx,vlist);
722 return (res);
723 case e_glassbox:
724 res = relman_glassbox_diffs(rel,filter,mtx);
725 return (res);
726 case e_opcode:
727 FPRINTF(stderr,"opcode differentiation not yet supported\n");
728 return 1.0e08;
729 case e_blackbox:
730 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
731 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
732 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
733 return res;
734 }
735 #endif
736 }
737
738 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
739 real64 relman_diffs_orig( struct rel_relation *rel, var_filter_t *filter,
740 mtx_matrix_t mtx)
741 {
742 real64 res = 0.0;
743 int32 row;
744 row = mtx_org_to_row(mtx,rel_sindex(rel));
745 if (rel_extnodeinfo(rel)) {
746 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
747 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
748 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
749 return res;
750 } else {
751 res -= exprman_diffs(rel,rel_rhs(rel),filter,row,mtx);
752 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
753 res += exprman_diffs(rel,rel_lhs(rel),filter,row,mtx);
754 return(res);
755 }
756 }
757 #endif
758
759 boolean relman_calc_satisfied( struct rel_relation *rel, real64 tolerance)
760 {
761 real64 res;
762 res = rel_residual(rel);
763 if (!finite(res)) {
764 rel_set_satisfied(rel,FALSE);
765 return( rel_satisfied(rel) );
766 }
767 if( res < 0.0 ) {
768 if( rel_less(rel) ) {
769 rel_set_satisfied(rel,TRUE);
770 return( rel_satisfied(rel) );
771 }
772 if( rel_greater(rel) ) {
773 rel_set_satisfied(rel,FALSE);
774 return( rel_satisfied(rel) );
775 }
776 rel_set_satisfied(rel,(-res <= tolerance ));
777 return( rel_satisfied(rel) );
778 }
779 if( res > 0.0 ) {
780 if( rel_greater(rel) ) {
781 rel_set_satisfied(rel,TRUE);
782 return( rel_satisfied(rel) );
783 }
784 if( rel_less(rel) ) {
785 rel_set_satisfied(rel,FALSE);
786 return( rel_satisfied(rel) );
787 }
788 rel_set_satisfied(rel,(res <= tolerance ));
789 return( rel_satisfied(rel) );
790 }
791 rel_set_satisfied(rel,rel_equal(rel)); /* strict >0 or <0 not satisfied */
792 return( rel_satisfied(rel) );
793 }
794
795 boolean relman_calc_satisfied_scaled(struct rel_relation *rel, real64 tolerance)
796 {
797 real64 res;
798 res = rel_residual(rel);
799
800 if( res < 0.0 )
801 if( rel_less(rel) )
802 rel_set_satisfied(rel,TRUE);
803 else if( rel_greater(rel) )
804 rel_set_satisfied(rel,FALSE);
805 else
806 rel_set_satisfied(rel,(-res <= tolerance * rel_nominal(rel)));
807 else if( res > 0.0 )
808 if( rel_greater(rel) )
809 rel_set_satisfied(rel,TRUE);
810 else if( rel_less(rel) )
811 rel_set_satisfied(rel,FALSE);
812 else
813 rel_set_satisfied(rel,(res <= tolerance * rel_nominal(rel)));
814 else
815 rel_set_satisfied(rel,rel_equal(rel));
816 return( rel_satisfied(rel) );
817 }
818
819 #if REIMPLEMENT
820 real64 *relman_directly_solve_new( struct rel_relation *rel,
821 struct var_variable *solvefor, int *able, int *nsolns,
822 real64 tolerance)
823 {
824 double *value;
825 if( rel_less(rel) || rel_greater(rel) || !rel_equal(rel) ||
826 rel_extnodeinfo(rel)) {
827 *able = FALSE;
828 *nsolns = 0;
829 return(NULL);
830 }
831 else if (rel->type == e_glassbox) {
832 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
833 return value;
834 }
835 else{
836 value =
837 exprman_directly_solve(rel,rel_lhs(rel),rel_rhs(rel),
838 solvefor,able,nsolns);
839 return value;
840 }
841 }
842
843 #else /* temporary */
844 real64 *relman_directly_solve_new( struct rel_relation *rel,
845 struct var_variable *solvefor, int *able, int *nsolns, real64 tolerance)
846 {
847 double *value;
848 if( rel_less(rel) ||
849 rel_greater(rel) ||
850 !rel_equal(rel) ||
851 rel_extnodeinfo(rel)) {
852 *able = FALSE;
853 *nsolns = 0;
854 return(NULL);
855 }
856 switch (rel->type ) {
857 case e_rel_glassbox:
858 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
859 return value;
860 case e_rel_token:
861 {
862 int nvars,n;
863 const struct var_variable **vlist;
864 unsigned long vindex; /* index to the compiler */
865 nvars = rel_n_incidences(rel);
866 vlist = rel_incidence_list(rel);
867 vindex = 0;
868 for (n=0;n <nvars; n++) {
869 if (vlist[n]==solvefor) {
870 vindex = n+1; /*compiler counts from 1 */
871 break;
872 }
873 }
874 value = RelationFindRoots(IPTR(rel_instance(rel)),
875 var_lower_bound(solvefor),
876 var_upper_bound(solvefor),
877 var_nominal(solvefor),
878 tolerance,
879 &(vindex),
880 able,nsolns);
881 return value;
882 }
883 default: /* e_rel_blackbox */
884 *able = FALSE;
885 *nsolns = 0;
886 return(NULL);
887 }
888 }
889 #endif
890
891 char *dummyrelstring(slv_system_t sys, struct rel_relation *rel, int style)
892 {
893 char *result;
894 UNUSED_PARAMETER(sys);
895 UNUSED_PARAMETER(rel);
896 UNUSED_PARAMETER(style);
897 result = ASC_NEW_ARRAY(char,80);
898 sprintf(result,"relman_make_*string_*fix not implemented yet");
899 return result;
900 }
901
902 /*
903 * converst the relations vlist into an array of int indices suitable
904 * for relation write string (compiler side indexing from 1, remember).
905 * if mORs == TRUE master indices are used, else solver.
906 */
907 static
908 void relman_cookup_indices(struct RXNameData *rd,
909 struct rel_relation *rel,int mORs)
910 {
911 int nvars,n;
912 const struct var_variable **vlist;
913
914 nvars = rel_n_incidences(rel);
915 rd->indices = rel_tmpalloc_array(1+nvars,int);
916 if (rd->indices == NULL) {
917 return;
918 }
919 vlist = rel_incidence_list(rel);
920 if (mORs) {
921 for (n = 0; n < nvars; n++) {
922 rd->indices[n+1] = var_mindex(vlist[n]);
923 }
924 } else {
925 for (n = 0; n < nvars; n++) {
926 rd->indices[n+1] = var_sindex(vlist[n]);
927 }
928 }
929 }
930
931 char *relman_make_vstring_infix(slv_system_t sys,
932 struct rel_relation *rel, int style)
933 {
934 char *sbeg;
935 struct RXNameData rd = {"x",NULL,""};
936
937 if (style) {
938 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
939 NULL,NULL,relio_ascend,NULL);
940 } else {
941 /* need to cook up output indices */
942 relman_cookup_indices(&rd,rel,1);
943 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
944 (WRSNameFunc)RelationVarXName,
945 &rd,relio_ascend,NULL);
946 }
947 return(sbeg);
948 }
949
950 char *relman_make_vstring_postfix(slv_system_t sys,
951 struct rel_relation *rel, int style)
952 {
953 char *sbeg;
954
955 if (style) {
956 sbeg = WriteRelationPostfixString(rel_instance(rel),slv_instance(sys));
957 } else {
958 #if REIMPLEMENT
959 left = exprman_make_xstring_postfix(rel,sys,rel_lhs(rel));
960 right = exprman_make_xstring_postfix(rel,sys,rel_rhs(rel));
961 #else
962 sbeg = ASC_NEW_ARRAY(char,60);
963 if (sbeg==NULL) return sbeg;
964 sprintf(sbeg,"relman_make_xstring_postfix not reimplemented.");
965 #endif
966 }
967 return(sbeg);
968 }

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