/[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 731 - (show annotations) (download) (as text)
Tue Jul 4 07:42:06 2006 UTC (14 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 26899 byte(s)
Removed some debug messages.
Fixed up return values for Integrators functions to comply with integrator.c API.
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 /*
393 CONSOLE_DEBUG("ABOUT TO CALL EXTREL_EVALUATE_RESIDUAL");
394 CONSOLE_DEBUG("REL_RELATION = %p",rel);
395 */
396 *calc_ok = 1;
397 res = ExtRel_Evaluate_Residual(rel);
398 rel_set_residual(rel,res);
399 return res;
400 }
401
402 /*
403 other types of relations (which include...). This latter approach is
404 apparently older and "reasonably correct, if slow".
405 */
406 /* CONSOLE_DEBUG("EVALUATE REL_RELATION = %p",rel); */
407 if(safe){
408 *calc_ok = (int32)RelationCalcResidualSafe(rel_instance(rel),&res);
409 safe_error_to_stderr( (enum safe_err *)calc_ok );
410
411 /* always set the relation residual when using safe functions */
412 rel_set_residual(rel,res);
413 }else{
414 *calc_ok = RelationCalcResidual(rel_instance(rel),&res);
415
416 /* for unsafe functions, set a fallback value in case of error */
417 if(*calc_ok){
418 res = 1.0e8;
419 }else{
420 rel_set_residual(rel,res);
421 }
422 }
423
424 /* flip the status flag: all values other than safe_ok become 0 */
425 *calc_ok = !(*calc_ok);
426 return res;
427 }
428
429 int32 relman_obj_direction(struct rel_relation *rel)
430 {
431 assert(rel!=NULL);
432 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
433 case e_minimize:
434 return -1;
435 case e_maximize:
436 return 1;
437 default:
438 return 0;
439 }
440 }
441
442 real64 relman_scale(struct rel_relation *rel)
443 {
444 real64 relnom;
445 assert(rel!=NULL);
446 relnom = CalcRelationNominal(rel_instance(rel));
447 if(relnom < 0.0000001) {
448 /* take care of small relnoms and relnom = 0 error returns */
449 relnom = 1.0;
450 }
451 rel_set_nominal(rel,relnom);
452 return relnom;
453 }
454
455 #if REIMPLEMENT /* compiler */
456 real64 relman_diff(struct rel_relation *rel, struct var_variable *var,
457 int safe)
458 {
459 /* FIX FIX FIX meaning kirk couldn't be botghered... */
460 real64 res = 0.0;
461 switch(rel->type) {
462 case e_glassbox:
463 case e_blackbox:
464 FPRINTF(stderr,"relman_diff not yet supported for black and glassbox\n");
465 return 1.0e08;
466 }
467 res -= exprman_diff(rel,rel_rhs(rel),var);
468 res += exprman_diff(rel,rel_lhs(rel),var);
469 return( res );
470 }
471 #endif
472
473 int relman_diff2(struct rel_relation *rel, var_filter_t *filter,
474 real64 *derivatives, int32 *variables,
475 int32 *count, int32 safe)
476 {
477 const struct var_variable **vlist=NULL;
478 real64 *gradient;
479 int32 len,c;
480 int status;
481
482 assert(rel!=NULL && filter!=NULL);
483 len = rel_n_incidences(rel);
484 vlist = rel_incidence_list(rel);
485
486 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
487 assert(gradient !=NULL);
488 *count = 0;
489 if( safe ) {
490 status =(int32)RelationCalcGradientSafe(rel_instance(rel),gradient);
491 safe_error_to_stderr( (enum safe_err *)&status );
492 /* always map when using safe functions */
493 for (c=0; c < len; c++) {
494 if (var_apply_filter(vlist[c],filter)) {
495 variables[*count] = var_sindex(vlist[c]);
496 derivatives[*count] = gradient[c];
497 (*count)++;
498 }
499 }
500 }
501 else {
502 if((status=RelationCalcGradient(rel_instance(rel),gradient)) == 0) {
503 /* successful */
504 for (c=0; c < len; c++) {
505 if (var_apply_filter(vlist[c],filter)) {
506 variables[*count] = var_sindex(vlist[c]);
507 derivatives[*count] = gradient[c];
508 (*count)++;
509 }
510 }
511 }
512 }
513
514 /* flip the status flag */
515 return !status;
516 }
517
518 int relman_diff_grad(struct rel_relation *rel, var_filter_t *filter,
519 real64 *derivatives, int32 *variables_master,
520 int32 *variables_solver, int32 *count, real64 *resid,
521 int32 safe)
522 {
523 const struct var_variable **vlist=NULL;
524 real64 *gradient;
525 int32 len,c;
526 int status;
527
528 assert(rel!=NULL && filter!=NULL);
529 len = rel_n_incidences(rel);
530 vlist = rel_incidence_list(rel);
531
532 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
533 assert(gradient !=NULL);
534 *count = 0;
535 if( safe ) {
536 /* CONSOLE_DEBUG("..."); */
537 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),
538 resid,gradient);
539 safe_error_to_stderr( (enum safe_err *)&status );
540 /* always map when using safe functions */
541 for (c=0; c < len; c++) {
542 if (var_apply_filter(vlist[c],filter)) {
543 variables_master[*count] = var_mindex(vlist[c]);
544 variables_solver[*count] = var_sindex(vlist[c]);
545 derivatives[*count] = gradient[c];
546 (*count)++;
547 }
548 }
549 }
550 else {
551 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient))== 0) {
552 /* successful */
553 for (c=0; c < len; c++) {
554 if (var_apply_filter(vlist[c],filter)) {
555 variables_master[*count] = var_mindex(vlist[c]);
556 variables_solver[*count] = var_sindex(vlist[c]);
557 derivatives[*count] = gradient[c];
558 (*count)++;
559 }
560 }
561 }
562 }
563
564 return !status; /* flip the status flag */
565 }
566
567 int32 relman_diff_harwell(struct rel_relation **rlist,
568 var_filter_t *vfilter, rel_filter_t *rfilter,
569 int32 rlen, int32 bias, int32 mORs,
570 real64 *avec, int32 *ivec, int32 *jvec)
571 {
572 const struct var_variable **vlist = NULL;
573 struct rel_relation *rel;
574 real64 residual, *resid;
575 real64 *gradient;
576 mtx_coord_t coord;
577 int32 len,c,r,k;
578 int32 errcnt;
579 enum safe_err status;
580
581 if (rlist == NULL || vfilter == NULL || rfilter == NULL || avec == NULL ||
582 rlen < 0 || mORs >3 || mORs < 0 || bias <0 || bias > 1) {
583 return 1;
584 }
585 resid = &residual;
586 errcnt = k = 0;
587
588 if ( ivec == NULL || jvec == NULL ) {
589 /*_skip stuffing ivec,jvec */
590 for (r=0; r < rlen; r++) {
591 rel = rlist[r];
592 len = rel_n_incidences(rel);
593 vlist = rel_incidence_list(rel);
594 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
595 if (gradient == NULL) {
596 return 1;
597 }
598 /* CONSOLE_DEBUG("..."); */
599 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
600 safe_error_to_stderr(&status);
601 if (status) {
602 errcnt--;
603 }
604 for (c=0; c < len; c++) {
605 if (var_apply_filter(vlist[c],vfilter)) {
606 avec[k] = gradient[c];
607 k++;
608 }
609 }
610 }
611 } else {
612 for (r=0; r < rlen; r++) {
613 rel = rlist[r];
614 len = rel_n_incidences(rel);
615 vlist = rel_incidence_list(rel);
616 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
617 if (gradient == NULL) {
618 return 1;
619 }
620 /* CONSOLE_DEBUG("..."); */
621 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
622 safe_error_to_stderr(&status);
623 if (status) {
624 errcnt--;
625 }
626 if (mORs & 2) {
627 coord.row = rel_sindex(rel);
628 } else {
629 coord.row = rel_mindex(rel);
630 }
631 for (c=0; c < len; c++) {
632 if (var_apply_filter(vlist[c],vfilter)) {
633 if (mORs & 1) {
634 coord.col = var_sindex(vlist[c]);
635 } else {
636 coord.col = var_mindex(vlist[c]);
637 }
638 avec[k] = gradient[c];
639 ivec[k] = coord.row;
640 jvec[k] = coord.col;
641 k++;
642 }
643 }
644 }
645 }
646 return errcnt;
647 }
648
649 int32 relman_jacobian_count(struct rel_relation **rlist, int32 rlen,
650 var_filter_t *vfilter,
651 rel_filter_t *rfilter, int32 *max)
652 {
653 int32 len, result=0, row, count, c;
654 const struct var_variable **list;
655 struct rel_relation *rel;
656 *max = 0;
657
658 for( row = 0; row < rlen; row++ ) {
659 rel = rlist[row];
660 if (rel_apply_filter(rel,rfilter)) {
661 len = rel_n_incidences(rel);
662 list = rel_incidence_list(rel);
663 count = 0;
664 for (c=0; c < len; c++) {
665 if( var_apply_filter(list[c],vfilter) ) {
666 count++;
667 }
668 }
669 result += count;
670 *max = MAX(*max,count);
671 }
672 }
673 return result;
674 }
675
676 int relman_diffs(struct rel_relation *rel, var_filter_t *filter,
677 mtx_matrix_t mtx, real64 *resid, int safe)
678 {
679 const struct var_variable **vlist=NULL;
680 real64 *gradient;
681 int32 len,c;
682 mtx_coord_t coord;
683 int status;
684
685 assert(rel!=NULL && filter!=NULL && mtx != NULL);
686 len = rel_n_incidences(rel);
687 vlist = rel_incidence_list(rel);
688 coord.row = rel_sindex(rel);
689 assert(coord.row>=0 && coord.row < mtx_order(mtx));
690
691 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
692 assert(gradient !=NULL);
693
694 /** @TODO fix this (it should all be in the compiler, or something) */
695 if(rel->nodeinfo){
696 /* CONSOLE_DEBUG("EVALUTING BLACKBOX DERIVATIVES FOR ROW %d",coord.row); */
697 *resid = extrel_resid_and_jacobian(rel, filter, coord.row, mtx);
698 return 0;
699 }
700
701 if(safe){
702 /* CONSOLE_DEBUG("..."); */
703 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
704 safe_error_to_stderr( (enum safe_err *)&status );
705 /* always map when using safe functions */
706 for (c=0; c < len; c++) {
707 if (var_apply_filter(vlist[c],filter)) {
708 coord.col = var_sindex(vlist[c]);
709 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
710 mtx_fill_org_value(mtx,&coord,gradient[c]);
711 }
712 }
713 }else{
714 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient)) == 0) {
715 /* successful */
716 for (c=0; c < len; c++) {
717 if (var_apply_filter(vlist[c],filter)) {
718 coord.col = var_sindex(vlist[c]);
719 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
720 mtx_fill_org_value(mtx,&coord,gradient[c]);
721 }
722 }
723 }
724 }
725
726 /* flip the status flag */
727 return !status;
728
729 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
730 switch (rel->type) {
731 case e_token:
732 vlist = rel_incidence_list(rel);
733 res -= exprman_diffs_alt(rel,rel_rhs(rel),filter,row,mtx,vlist);
734 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
735 res += exprman_diffs_alt(rel,rel_lhs(rel),filter,row,mtx,vlist);
736 return (res);
737 case e_glassbox:
738 res = relman_glassbox_diffs(rel,filter,mtx);
739 return (res);
740 case e_opcode:
741 FPRINTF(stderr,"opcode differentiation not yet supported\n");
742 return 1.0e08;
743 case e_blackbox:
744 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
745 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
746 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
747 return res;
748 }
749 #endif
750 }
751
752 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
753 real64 relman_diffs_orig( struct rel_relation *rel, var_filter_t *filter,
754 mtx_matrix_t mtx)
755 {
756 real64 res = 0.0;
757 int32 row;
758 row = mtx_org_to_row(mtx,rel_sindex(rel));
759 if (rel_extnodeinfo(rel)) {
760 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
761 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
762 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
763 return res;
764 } else {
765 res -= exprman_diffs(rel,rel_rhs(rel),filter,row,mtx);
766 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
767 res += exprman_diffs(rel,rel_lhs(rel),filter,row,mtx);
768 return(res);
769 }
770 }
771 #endif
772
773 boolean relman_calc_satisfied( struct rel_relation *rel, real64 tolerance)
774 {
775 real64 res;
776 res = rel_residual(rel);
777 if (!finite(res)) {
778 rel_set_satisfied(rel,FALSE);
779 return( rel_satisfied(rel) );
780 }
781 if( res < 0.0 ) {
782 if( rel_less(rel) ) {
783 rel_set_satisfied(rel,TRUE);
784 return( rel_satisfied(rel) );
785 }
786 if( rel_greater(rel) ) {
787 rel_set_satisfied(rel,FALSE);
788 return( rel_satisfied(rel) );
789 }
790 rel_set_satisfied(rel,(-res <= tolerance ));
791 return( rel_satisfied(rel) );
792 }
793 if( res > 0.0 ) {
794 if( rel_greater(rel) ) {
795 rel_set_satisfied(rel,TRUE);
796 return( rel_satisfied(rel) );
797 }
798 if( rel_less(rel) ) {
799 rel_set_satisfied(rel,FALSE);
800 return( rel_satisfied(rel) );
801 }
802 rel_set_satisfied(rel,(res <= tolerance ));
803 return( rel_satisfied(rel) );
804 }
805 rel_set_satisfied(rel,rel_equal(rel)); /* strict >0 or <0 not satisfied */
806 return( rel_satisfied(rel) );
807 }
808
809 boolean relman_calc_satisfied_scaled(struct rel_relation *rel, real64 tolerance)
810 {
811 real64 res;
812 res = rel_residual(rel);
813
814 if( res < 0.0 )
815 if( rel_less(rel) )
816 rel_set_satisfied(rel,TRUE);
817 else if( rel_greater(rel) )
818 rel_set_satisfied(rel,FALSE);
819 else
820 rel_set_satisfied(rel,(-res <= tolerance * rel_nominal(rel)));
821 else if( res > 0.0 )
822 if( rel_greater(rel) )
823 rel_set_satisfied(rel,TRUE);
824 else if( rel_less(rel) )
825 rel_set_satisfied(rel,FALSE);
826 else
827 rel_set_satisfied(rel,(res <= tolerance * rel_nominal(rel)));
828 else
829 rel_set_satisfied(rel,rel_equal(rel));
830 return( rel_satisfied(rel) );
831 }
832
833 #if REIMPLEMENT
834 real64 *relman_directly_solve_new( struct rel_relation *rel,
835 struct var_variable *solvefor, int *able, int *nsolns,
836 real64 tolerance)
837 {
838 double *value;
839 if( rel_less(rel) || rel_greater(rel) || !rel_equal(rel) ||
840 rel_extnodeinfo(rel)) {
841 *able = FALSE;
842 *nsolns = 0;
843 return(NULL);
844 }
845 else if (rel->type == e_glassbox) {
846 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
847 return value;
848 }
849 else{
850 value =
851 exprman_directly_solve(rel,rel_lhs(rel),rel_rhs(rel),
852 solvefor,able,nsolns);
853 return value;
854 }
855 }
856
857 #else /* temporary */
858 real64 *relman_directly_solve_new( struct rel_relation *rel,
859 struct var_variable *solvefor, int *able, int *nsolns, real64 tolerance)
860 {
861 double *value;
862 if( rel_less(rel) ||
863 rel_greater(rel) ||
864 !rel_equal(rel) ||
865 rel_extnodeinfo(rel)) {
866 *able = FALSE;
867 *nsolns = 0;
868 return(NULL);
869 }
870 switch (rel->type ) {
871 case e_rel_glassbox:
872 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
873 return value;
874 case e_rel_token:
875 {
876 int nvars,n;
877 const struct var_variable **vlist;
878 unsigned long vindex; /* index to the compiler */
879 nvars = rel_n_incidences(rel);
880 vlist = rel_incidence_list(rel);
881 vindex = 0;
882 for (n=0;n <nvars; n++) {
883 if (vlist[n]==solvefor) {
884 vindex = n+1; /*compiler counts from 1 */
885 break;
886 }
887 }
888 value = RelationFindRoots(IPTR(rel_instance(rel)),
889 var_lower_bound(solvefor),
890 var_upper_bound(solvefor),
891 var_nominal(solvefor),
892 tolerance,
893 &(vindex),
894 able,nsolns);
895 return value;
896 }
897 default: /* e_rel_blackbox */
898 *able = FALSE;
899 *nsolns = 0;
900 return(NULL);
901 }
902 }
903 #endif
904
905 char *dummyrelstring(slv_system_t sys, struct rel_relation *rel, int style)
906 {
907 char *result;
908 UNUSED_PARAMETER(sys);
909 UNUSED_PARAMETER(rel);
910 UNUSED_PARAMETER(style);
911 result = ASC_NEW_ARRAY(char,80);
912 sprintf(result,"relman_make_*string_*fix not implemented yet");
913 return result;
914 }
915
916 /*
917 * converst the relations vlist into an array of int indices suitable
918 * for relation write string (compiler side indexing from 1, remember).
919 * if mORs == TRUE master indices are used, else solver.
920 */
921 static
922 void relman_cookup_indices(struct RXNameData *rd,
923 struct rel_relation *rel,int mORs)
924 {
925 int nvars,n;
926 const struct var_variable **vlist;
927
928 nvars = rel_n_incidences(rel);
929 rd->indices = rel_tmpalloc_array(1+nvars,int);
930 if (rd->indices == NULL) {
931 return;
932 }
933 vlist = rel_incidence_list(rel);
934 if (mORs) {
935 for (n = 0; n < nvars; n++) {
936 rd->indices[n+1] = var_mindex(vlist[n]);
937 }
938 } else {
939 for (n = 0; n < nvars; n++) {
940 rd->indices[n+1] = var_sindex(vlist[n]);
941 }
942 }
943 }
944
945 char *relman_make_vstring_infix(slv_system_t sys,
946 struct rel_relation *rel, int style)
947 {
948 char *sbeg;
949 struct RXNameData rd = {"x",NULL,""};
950
951 if (style) {
952 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
953 NULL,NULL,relio_ascend,NULL);
954 } else {
955 /* need to cook up output indices */
956 relman_cookup_indices(&rd,rel,1);
957 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
958 (WRSNameFunc)RelationVarXName,
959 &rd,relio_ascend,NULL);
960 }
961 return(sbeg);
962 }
963
964 char *relman_make_vstring_postfix(slv_system_t sys,
965 struct rel_relation *rel, int style)
966 {
967 char *sbeg;
968
969 if (style) {
970 sbeg = WriteRelationPostfixString(rel_instance(rel),slv_instance(sys));
971 } else {
972 #if REIMPLEMENT
973 left = exprman_make_xstring_postfix(rel,sys,rel_lhs(rel));
974 right = exprman_make_xstring_postfix(rel,sys,rel_rhs(rel));
975 #else
976 sbeg = ASC_NEW_ARRAY(char,60);
977 if (sbeg==NULL) return sbeg;
978 sprintf(sbeg,"relman_make_xstring_postfix not reimplemented.");
979 #endif
980 }
981 return(sbeg);
982 }

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