/[ascend]/trunk/ascend4/solver/relman.c
ViewVC logotype

Contents of /trunk/ascend4/solver/relman.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 26795 byte(s)
Setting up web subdirectory in repository
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 "compiler/extfunc.h"
41 #include "compiler/dimen.h"
42 #include "compiler/types.h"
43 #include "compiler/find.h"
44 #include "general/list.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 "solver/mtx.h"
53 #include "solver/slv_types.h"
54 #include "solver/var.h"
55 #include "solver/rel.h"
56 #include "solver/discrete.h"
57 #include "solver/conditional.h"
58 #include "solver/logrel.h"
59 #include "solver/bnd.h"
60 #include "solver/relman.h"
61 #include "compiler/rootfind.h"
62 #include "solver/slv_server.h"
63
64 #define IPTR(i) ((struct Instance *)(i))
65
66 #define KILL 0 /* compile dead code if kill = 1 */
67 #define REIMPLEMENT 0 /* code that needs to be reimplemented */
68
69 static POINTER rel_tmpalloc( int nbytes)
70 /**
71 *** Temporarily allocates a given number of bytes. The memory need
72 *** not be freed, but the next call to this function will reuse the
73 *** previous allocation. Memory returned will NOT be zeroed.
74 *** Calling with nbytes==0 will free any memory allocated.
75 **/
76 {
77 static char *ptr = NULL;
78 static int cap = 0;
79
80 if (nbytes) {
81 if( nbytes > cap ) {
82 if( ptr != NULL ) ascfree(ptr);
83 ptr = (char *)ascmalloc(nbytes);
84 cap = nbytes;
85 }
86 } else {
87 if (ptr) ascfree(ptr);
88 ptr=NULL;
89 cap=0;
90 }
91 if ( cap >0) {
92 return(ptr);
93 } else {
94 return NULL;
95 }
96 }
97
98 #define rel_tmpalloc_array(nelts,type) \
99 ((nelts) > 0 ? (type *)tmpalloc((nelts)*sizeof(type)) : NULL)
100 /**
101 *** Creates an array of "nelts" objects, each with type "type".
102 **/
103
104 void relman_free_reused_mem(void)
105 {
106 /* rel_tmpalloc(0); */
107 RelationFindRoots(NULL,0,0,0,0,NULL,NULL,NULL);
108 }
109
110
111 #if REIMPLEMENT
112 boolean relman_is_linear( struct rel_relation *rel, var_filter_t *filter)
113 {
114 return (
115 exprman_is_linear(rel,rel_lhs(rel),filter) &&
116 exprman_is_linear(rel,rel_rhs(rel),filter)
117 );
118 }
119
120 real64 relman_linear_coef(struct rel_relation *rel, struct var_variable *var,
121 var_filter_t *filter)
122 {
123 return(
124 exprman_linear_coef(rel,rel_lhs(rel),var,filter) -
125 exprman_linear_coef(rel,rel_rhs(rel),var,filter)
126 );
127 }
128 #endif
129
130 #if KILL
131 void relman_decide_incidence( struct rel_relation *rel)
132 {
133 struct var_variable **list;
134 int c;
135
136 list = rel_incidence_list(rel);
137 for( c=0; list[c] != NULL; c++ )
138 var_set_incident(list[c],TRUE);
139 }
140 #endif
141
142 void relman_get_incidence(struct rel_relation *rel, var_filter_t *filter,
143 mtx_matrix_t mtx)
144 {
145 const struct var_variable **list;
146 mtx_coord_t nz;
147 int c,len;
148
149 assert(rel!=NULL && filter !=NULL && mtx != NULL);
150 nz.row = rel_sindex(rel);
151 len = rel_n_incidences(rel);
152
153 list = rel_incidence_list(rel);
154 for (c=0; c < len; c++) {
155 if( var_apply_filter(list[c],filter) ) {
156 nz.col = var_sindex(list[c]);
157 mtx_fill_org_value(mtx,&nz,1.0);
158 }
159 }
160 }
161
162 /*
163 *********************************************************************
164 * Code to deal with glassbox relations processing.
165 *********************************************************************
166 */
167
168 static double dsolve_scratch = 0.0; /* some workspace */
169 #define DSOLVE_TOLERANCE 1.0e-08 /* no longer needed */
170
171 static
172 real64 *relman_glassbox_dsolve(struct rel_relation *rel,
173 struct var_variable *solvefor,
174 int *able,
175 int *nsolns,
176 real64 tolerance)
177 {
178 int n,m,j,dummy = 0;
179 int mode, result;
180 int index,solve_for_index = -1;
181 struct Instance *var;
182 CONST struct relation *cmplr_reln;
183 enum Expr_enum type;
184 CONST struct gl_list_t *vlist;
185 double *f, *x, *g, value;
186 double lower, upper, nominal;
187 ExtEvalFunc **table;
188
189 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
190 assert(type==e_glassbox);
191
192 vlist = RelationVarList(cmplr_reln);
193 m = rel_sindex(rel);
194 n = (int)gl_length(vlist);
195 f = (real64 *)asccalloc((1 + 2*n),sizeof(double));
196 x = &f[1];
197 g = &f[n+1];
198
199 /*
200 * Load x vector, and get the index into the x[vector]
201 * of the variable that we are solving for.
202 */
203 for (j=0;j<n;j++) {
204 var = (struct Instance *)gl_fetch(vlist,(unsigned long)(j+1));
205 x[j] = RealAtomValue(var);
206 if (var_instance(solvefor)== var)
207 solve_for_index = j;
208 }
209
210 /*
211 * Set up bounds and tolerance.
212 */
213 lower = var_lower_bound(solvefor);
214 upper = var_upper_bound(solvefor);
215 nominal = var_nominal(solvefor);
216 /*tolerance = DSOLVE_TOLERANCE; now passed in. */
217
218 /*
219 * Get the evaluation function.
220 */
221 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
222 index = GlassBoxRelIndex(cmplr_reln);
223
224 /*
225 * Call the rootfinder. A note of CAUTION:
226 * zbrent is going to call the evaluation function.
227 * It expects, that the results of this evaluation will be
228 * written to f[m]. GlassBox relations however always
229 * write to slot f[0]. To keep the world happy we send in
230 * a dummy = 0. This needs to be made cleaner. !! FIX !!
231 */
232 value = zbrent(table[index], &lower, &upper,
233 &mode, &dummy, &solve_for_index,
234 x, x , f, g, &tolerance, &result);
235 if (result==0) {
236 dsolve_scratch = value;
237 ascfree((char *)f);
238 *able = TRUE;
239 *nsolns = 1;
240 return &dsolve_scratch;
241 }
242 else{
243 dsolve_scratch = 0.0;
244 ascfree((char *)f);
245 *able = FALSE;
246 *nsolns = 0;
247 return NULL;
248 }
249 }
250
251
252 #ifdef THIS_IS_AN_UNUSED_FUNCTION
253 static
254 real64 relman_glassbox_eval(struct rel_relation *rel)
255 {
256 int n,m,mode,result;
257 int index,j;
258 CONST struct relation *cmplr_reln;
259 enum Expr_enum type;
260 CONST struct gl_list_t *vlist;
261 double *f, *x, *g, value;
262 ExtEvalFunc **table, *eval;
263
264 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
265 assert(type==e_glassbox);
266
267 vlist = RelationVarList(cmplr_reln);
268 m = rel_sindex(rel);
269 n = (int)gl_length(vlist);
270 /* this needs to be reused memory, fergossake */
271 f = (real64 *)asccalloc((1 + 2*n),sizeof(double)); /* resid */
272 x = &f[1]; /* var values */
273 g = &f[n+1]; /* gradient */
274
275 for (j=0;j<n;j++) {
276 x[j] = RealAtomValue(gl_fetch(vlist,(unsigned long)(j+1)));
277 }
278
279 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
280 index = GlassBoxRelIndex(cmplr_reln);
281
282 /*
283 * Remember to set up the fpe traps etc ....
284 * and to monitor the return flag.
285 */
286 mode = 0;
287 eval = table[index];
288 result = (*eval)(&mode,&m,&n,x,NULL,f,g);
289 value = f[0];
290 ascfree((char *)f);
291 return value;
292 }
293 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
294
295
296 #define BROKENKIRK 0
297 #if BROKENKIRK
298 /* fills filter passing gradient elements to matrix */
299 /* this needs to be buried on the compiler side. */
300 void relman_map_grad2mtx( struct rel_relation *rel, var_filter_t *filter,
301 mtx_matrix_t mtx, CONST struct gl_list_t *varlist, double *g, int m, int n)
302 {
303 mtx_coord_t nz;
304 struct var_variable *var;
305 double value;
306 int j;
307
308 nz.row = m; /* org row of glassbox reln? rel_sindex? */
309
310 for (j=0;j<n;j++) {
311 /* var = (struct Instance *)gl_fetch(varlist,(unsigned long)(j+1)); */
312 var = rel_incidence(rel)[j];
313 if (var_apply_filter(var)) {
314 nz.col = var_sindex(var);
315 value = g[j] /* + mtx_value(mtx,&nz) no longer incrementing row */;
316 mtx_fill_org_value(mtx,&nz, value);
317 }
318 }
319 }
320
321 real64 relman_glassbox_diffs( struct rel_relation *rel,
322 var_filter_t *filter, mtx_matrix_t mtx)
323 {
324 int n,m,mode,result;
325 int index,j;
326 struct Instance *var;
327 CONST struct relation *cmplr_reln;
328 enum Expr_enum type;
329 CONST struct gl_list_t *vlist;
330 double *f, *x, *g, value;
331 ExtEvalFunc **table, *eval;
332
333 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
334 vlist = RelationVarList(cmplr_reln);
335 m = rel_sindex(rel);
336 n = (int)gl_length(vlist);
337 /* this needs to be reused memory ! */
338 f = (real64 *)asccalloc((1 + 2*n),sizeof(double));
339 x = &f[1];
340 g = &f[n+1];
341
342 for (j=0;j<n;j++) {
343 x[j] = RealAtomValue(gl_fetch(vlist,(unsigned long)j+1));
344 }
345
346 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
347 index = GlassBoxRelIndex(cmplr_reln);
348 eval = table[index];
349 result = (*eval)(0,&m,&n,x,NULL,f,g);
350
351 table = GetDerivJumpTable(GlassBoxExtFunc(cmplr_reln));
352
353 mode = 0;
354 eval = table[index];
355 result += (*eval)(&mode,&m,&n,x,NULL,f,g);
356 relman_map_grad2mtx(rel,filter,mtx,vlist,g,m,n);
357
358 /*
359 * Remember to set up the fpe traps etc ....
360 * and to monitor the return flag.
361 */
362 value = f[0];
363 ascfree((char *)f);
364 return value;
365 }
366
367 #else
368 #define relman_glassbox_diffs(rel,filter,mtx) abort()
369 #endif
370
371 real64 relman_eval(struct rel_relation *rel, int32 *status, int safe)
372 {
373 real64 res;
374 assert(status!=NULL && rel!=NULL);
375 if ( rel->type == e_rel_token ) {
376 if (!RelationCalcResidualBinary(
377 GetInstanceRelationOnly(IPTR(rel->instance)),&res)) {
378 *status = 1; /* calc_ok */
379 rel_set_residual(rel,res);
380 return res;
381 }
382 /* else we don't care -- go on to the old handling which
383 * is reasonably correct, if slow.
384 */
385 }
386 if( safe ) {
387 *status = (int32)RelationCalcResidualSafe(rel_instance(rel),&res);
388 safe_error_to_stderr( (enum safe_err *)status );
389 /* always set the relation residual when using safe functions */
390 rel_set_residual(rel,res);
391 } else {
392 *status = RelationCalcResidual(rel_instance(rel),&res);
393 if ( *status ) {
394 /* an error occured */
395 res = 1.0e8;
396 } else {
397 rel_set_residual(rel,res);
398 }
399 }
400 /* flip the status flag: all values other than safe_ok become 0 */
401 *status = !(*status);
402 return res;
403
404 #if REIMPLEMENT /* all this needs to be done on the compiler side */
405 it may take some changes in the CalcResidual header to do so.
406 switch (rel->type) {
407 case e_token:
408 res = exprman_eval(rel,rel_lhs(rel)) - exprman_eval(rel,rel_rhs(rel));
409 break;
410 case e_opcode:
411 FPRINTF(stderr,"opcode relation processing not yet supported\n");
412 res = 1.0e08;
413 break;
414 case e_glassbox:
415 res = relman_glassbox_eval(rel);
416 break;
417 case e_blackbox:
418 res = ExtRel_Evaluate_LHS(rel) - ExtRel_Evaluate_RHS(rel);
419 break;
420 default:
421 FPRINTF(stderr,"unknown relation type in (relman_eval)\n");
422 res = 1.0e08;
423 break;
424 }
425 #endif
426
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 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),
537 resid,gradient);
538 safe_error_to_stderr( (enum safe_err *)&status );
539 /* always map when using safe functions */
540 for (c=0; c < len; c++) {
541 if (var_apply_filter(vlist[c],filter)) {
542 variables_master[*count] = var_mindex(vlist[c]);
543 variables_solver[*count] = var_sindex(vlist[c]);
544 derivatives[*count] = gradient[c];
545 (*count)++;
546 }
547 }
548 }
549 else {
550 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient))== 0) {
551 /* successful */
552 for (c=0; c < len; c++) {
553 if (var_apply_filter(vlist[c],filter)) {
554 variables_master[*count] = var_mindex(vlist[c]);
555 variables_solver[*count] = var_sindex(vlist[c]);
556 derivatives[*count] = gradient[c];
557 (*count)++;
558 }
559 }
560 }
561 }
562
563 return !status; /* flip the status flag */
564 }
565
566 int32 relman_diff_harwell(struct rel_relation **rlist,
567 var_filter_t *vfilter, rel_filter_t *rfilter,
568 int32 rlen, int32 bias, int32 mORs,
569 real64 *avec, int32 *ivec, int32 *jvec)
570 {
571 const struct var_variable **vlist = NULL;
572 struct rel_relation *rel;
573 real64 residual, *resid;
574 real64 *gradient;
575 mtx_coord_t coord;
576 int32 len,c,r,k;
577 int32 errcnt;
578 enum safe_err status;
579
580 if (rlist == NULL || vfilter == NULL || rfilter == NULL || avec == NULL ||
581 rlen < 0 || mORs >3 || mORs < 0 || bias <0 || bias > 1) {
582 return 1;
583 }
584 resid = &residual;
585 errcnt = k = 0;
586
587 if ( ivec == NULL || jvec == NULL ) {
588 /*_skip stuffing ivec,jvec */
589 for (r=0; r < rlen; r++) {
590 rel = rlist[r];
591 len = rel_n_incidences(rel);
592 vlist = rel_incidence_list(rel);
593 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
594 if (gradient == NULL) {
595 return 1;
596 }
597 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
598 safe_error_to_stderr(&status);
599 if (status) {
600 errcnt--;
601 }
602 for (c=0; c < len; c++) {
603 if (var_apply_filter(vlist[c],vfilter)) {
604 avec[k] = gradient[c];
605 k++;
606 }
607 }
608 }
609 } else {
610 for (r=0; r < rlen; r++) {
611 rel = rlist[r];
612 len = rel_n_incidences(rel);
613 vlist = rel_incidence_list(rel);
614 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
615 if (gradient == NULL) {
616 return 1;
617 }
618 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
619 safe_error_to_stderr(&status);
620 if (status) {
621 errcnt--;
622 }
623 if (mORs & 2) {
624 coord.row = rel_sindex(rel);
625 } else {
626 coord.row = rel_mindex(rel);
627 }
628 for (c=0; c < len; c++) {
629 if (var_apply_filter(vlist[c],vfilter)) {
630 if (mORs & 1) {
631 coord.col = var_sindex(vlist[c]);
632 } else {
633 coord.col = var_mindex(vlist[c]);
634 }
635 avec[k] = gradient[c];
636 ivec[k] = coord.row;
637 jvec[k] = coord.col;
638 k++;
639 }
640 }
641 }
642 }
643 return errcnt;
644 }
645
646 int32 relman_jacobian_count(struct rel_relation **rlist, int32 rlen,
647 var_filter_t *vfilter,
648 rel_filter_t *rfilter, int32 *max)
649 {
650 int32 len, result=0, row, count, c;
651 const struct var_variable **list;
652 struct rel_relation *rel;
653 *max = 0;
654
655 for( row = 0; row < rlen; row++ ) {
656 rel = rlist[row];
657 if (rel_apply_filter(rel,rfilter)) {
658 len = rel_n_incidences(rel);
659 list = rel_incidence_list(rel);
660 count = 0;
661 for (c=0; c < len; c++) {
662 if( var_apply_filter(list[c],vfilter) ) {
663 count++;
664 }
665 }
666 result += count;
667 *max = MAX(*max,count);
668 }
669 }
670 return result;
671 }
672
673 int relman_diffs(struct rel_relation *rel, var_filter_t *filter,
674 mtx_matrix_t mtx, real64 *resid, int safe)
675 {
676 const struct var_variable **vlist=NULL;
677 real64 *gradient;
678 int32 len,c;
679 mtx_coord_t coord;
680 int status;
681
682 assert(rel!=NULL && filter!=NULL && mtx != NULL);
683 len = rel_n_incidences(rel);
684 vlist = rel_incidence_list(rel);
685 coord.row = rel_sindex(rel);
686 assert(coord.row>=0 && coord.row < mtx_order(mtx));
687
688 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
689 assert(gradient !=NULL);
690 if( safe ) {
691 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
692 safe_error_to_stderr( (enum safe_err *)&status );
693 /* always map when using safe functions */
694 for (c=0; c < len; c++) {
695 if (var_apply_filter(vlist[c],filter)) {
696 coord.col = var_sindex(vlist[c]);
697 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
698 mtx_fill_org_value(mtx,&coord,gradient[c]);
699 }
700 }
701 }
702 else {
703 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient)) == 0) {
704 /* successful */
705 for (c=0; c < len; c++) {
706 if (var_apply_filter(vlist[c],filter)) {
707 coord.col = var_sindex(vlist[c]);
708 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
709 mtx_fill_org_value(mtx,&coord,gradient[c]);
710 }
711 }
712 }
713 }
714
715 /* flip the status flag */
716 return !status;
717
718 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
719 switch (rel->type) {
720 case e_token:
721 vlist = rel_incidence_list(rel);
722 res -= exprman_diffs_alt(rel,rel_rhs(rel),filter,row,mtx,vlist);
723 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
724 res += exprman_diffs_alt(rel,rel_lhs(rel),filter,row,mtx,vlist);
725 return (res);
726 case e_glassbox:
727 res = relman_glassbox_diffs(rel,filter,mtx);
728 return (res);
729 case e_opcode:
730 FPRINTF(stderr,"opcode differentiation not yet supported\n");
731 return 1.0e08;
732 case e_blackbox:
733 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
734 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
735 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
736 return res;
737 }
738 #endif
739 }
740
741 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
742 real64 relman_diffs_orig( struct rel_relation *rel, var_filter_t *filter,
743 mtx_matrix_t mtx)
744 {
745 real64 res = 0.0;
746 int32 row;
747 row = mtx_org_to_row(mtx,rel_sindex(rel));
748 if (rel_extnodeinfo(rel)) {
749 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
750 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
751 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
752 return res;
753 } else {
754 res -= exprman_diffs(rel,rel_rhs(rel),filter,row,mtx);
755 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
756 res += exprman_diffs(rel,rel_lhs(rel),filter,row,mtx);
757 return(res);
758 }
759 }
760 #endif
761
762 boolean relman_calc_satisfied( struct rel_relation *rel, real64 tolerance)
763 {
764 real64 res;
765 res = rel_residual(rel);
766 if (!finite(res)) {
767 rel_set_satisfied(rel,FALSE);
768 return( rel_satisfied(rel) );
769 }
770 if( res < 0.0 ) {
771 if( rel_less(rel) ) {
772 rel_set_satisfied(rel,TRUE);
773 return( rel_satisfied(rel) );
774 }
775 if( rel_greater(rel) ) {
776 rel_set_satisfied(rel,FALSE);
777 return( rel_satisfied(rel) );
778 }
779 rel_set_satisfied(rel,(-res <= tolerance ));
780 return( rel_satisfied(rel) );
781 }
782 if( res > 0.0 ) {
783 if( rel_greater(rel) ) {
784 rel_set_satisfied(rel,TRUE);
785 return( rel_satisfied(rel) );
786 }
787 if( rel_less(rel) ) {
788 rel_set_satisfied(rel,FALSE);
789 return( rel_satisfied(rel) );
790 }
791 rel_set_satisfied(rel,(res <= tolerance ));
792 return( rel_satisfied(rel) );
793 }
794 rel_set_satisfied(rel,rel_equal(rel)); /* strict >0 or <0 not satisfied */
795 return( rel_satisfied(rel) );
796 }
797
798 boolean relman_calc_satisfied_scaled(struct rel_relation *rel, real64 tolerance)
799 {
800 real64 res;
801 res = rel_residual(rel);
802
803 if( res < 0.0 )
804 if( rel_less(rel) )
805 rel_set_satisfied(rel,TRUE);
806 else if( rel_greater(rel) )
807 rel_set_satisfied(rel,FALSE);
808 else
809 rel_set_satisfied(rel,(-res <= tolerance * rel_nominal(rel)));
810 else if( res > 0.0 )
811 if( rel_greater(rel) )
812 rel_set_satisfied(rel,TRUE);
813 else if( rel_less(rel) )
814 rel_set_satisfied(rel,FALSE);
815 else
816 rel_set_satisfied(rel,(res <= tolerance * rel_nominal(rel)));
817 else
818 rel_set_satisfied(rel,rel_equal(rel));
819 return( rel_satisfied(rel) );
820 }
821
822 #if REIMPLEMENT
823 real64 *relman_directly_solve_new( struct rel_relation *rel,
824 struct var_variable *solvefor, int *able, int *nsolns,
825 real64 tolerance)
826 {
827 double *value;
828 if( rel_less(rel) || rel_greater(rel) || !rel_equal(rel) ||
829 rel_extnodeinfo(rel)) {
830 *able = FALSE;
831 *nsolns = 0;
832 return(NULL);
833 }
834 else if (rel->type == e_glassbox) {
835 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
836 return value;
837 }
838 else{
839 value =
840 exprman_directly_solve(rel,rel_lhs(rel),rel_rhs(rel),
841 solvefor,able,nsolns);
842 return value;
843 }
844 }
845
846 #else /* temporary */
847 real64 *relman_directly_solve_new( struct rel_relation *rel,
848 struct var_variable *solvefor, int *able, int *nsolns, real64 tolerance)
849 {
850 double *value;
851 if( rel_less(rel) ||
852 rel_greater(rel) ||
853 !rel_equal(rel) ||
854 rel_extnodeinfo(rel)) {
855 *able = FALSE;
856 *nsolns = 0;
857 return(NULL);
858 }
859 switch (rel->type ) {
860 case e_rel_glassbox:
861 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
862 return value;
863 case e_rel_token:
864 {
865 int nvars,n;
866 const struct var_variable **vlist;
867 unsigned long vindex; /* index to the compiler */
868 nvars = rel_n_incidences(rel);
869 vlist = rel_incidence_list(rel);
870 vindex = 0;
871 for (n=0;n <nvars; n++) {
872 if (vlist[n]==solvefor) {
873 vindex = n+1; /*compiler counts from 1 */
874 break;
875 }
876 }
877 value = RelationFindRoots(IPTR(rel_instance(rel)),
878 var_lower_bound(solvefor),
879 var_upper_bound(solvefor),
880 var_nominal(solvefor),
881 tolerance,
882 &(vindex),
883 able,nsolns);
884 return value;
885 }
886 default: /* e_rel_blackbox */
887 *able = FALSE;
888 *nsolns = 0;
889 return(NULL);
890 }
891 }
892 #endif
893
894 char *dummyrelstring(slv_system_t sys, struct rel_relation *rel, int style)
895 {
896 char *result;
897 result = (char *)ascmalloc(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 = (char *)ascmalloc(60);
963 if (sbeg==NULL) return sbeg;
964 sprintf(sbeg,"relman_make_xstring_postfix not reimplemented.");
965 #endif
966 }
967 return(sbeg);
968 }

Properties

Name Value
svn:executable *

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