/[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 1034 - (show annotations) (download) (as text)
Thu Jan 4 05:37:55 2007 UTC (13 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 26576 byte(s)
Switch 'void slv_*' function to output a 0-on-success error flag.
Added some debug output to ensure that such output gets through.
Fixed TestBlackBox.testfail1 to work, using above changes.
Simulation::solve now throws an exception on failure (will need to modify PyGTK GUI accordingly)
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3 Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
4 Copyright (C) 1993 Joseph Zaher
5 Copyright (C) 1990 Karl Michael Westerberg
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2, or (at your option)
10 any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 59 Temple Place - Suite 330,
20 Boston, MA 02111-1307, USA.
21 *//** @file
22 Relation manipulator module for the SLV solver.
23 *//*
24 by Karl Michael Westerberg. Created 2/6/90
25 Last in CVS: $Revision: 1.44 $ $Date: 1998/04/23 23:56:22 $ $Author: ballan $
26 */
27
28 #include <math.h>
29 #include <utilities/ascConfig.h>
30 #include <compiler/instance_enum.h>
31 #include <compiler/fractions.h>
32 #include <compiler/compiler.h>
33 #include <utilities/ascMalloc.h>
34 #include <compiler/functype.h>
35 #include <compiler/safe.h>
36 #include <general/list.h>
37 #include <compiler/extfunc.h>
38 #include <compiler/dimen.h>
39 #include <compiler/expr_types.h>
40 #include <compiler/find.h>
41 #include <compiler/atomvalue.h>
42 #include <compiler/mathinst.h>
43 #include <compiler/relation_type.h>
44 #include <compiler/rel_blackbox.h>
45 #include <compiler/vlist.h>
46 #include <compiler/relation.h>
47 #include <compiler/relation_util.h>
48 #include <compiler/relation_io.h>
49 #define _SLV_SERVER_C_SEEN_
50 #include "mtx.h"
51 #include "slv_types.h"
52 #include "var.h"
53 #include "rel.h"
54 #include "discrete.h"
55 #include "conditional.h"
56 #include "logrel.h"
57 #include "bnd.h"
58 #include "relman.h"
59 #include <compiler/rootfind.h>
60 #include "slv_server.h"
61 #include <general/mathmacros.h>
62
63 #define IPTR(i) ((struct Instance *)(i))
64
65 #define KILL 0 /* compile dead code if kill = 1 */
66 #define REIMPLEMENT 0 /* code that needs to be reimplemented */
67
68 static POINTER rel_tmpalloc( int nbytes)
69 /**
70 *** Temporarily allocates a given number of bytes. The memory need
71 *** not be freed, but the next call to this function will reuse the
72 *** previous allocation. Memory returned will NOT be zeroed.
73 *** Calling with nbytes==0 will free any memory allocated.
74 **/
75 {
76 static char *ptr = NULL;
77 static int cap = 0;
78
79 if (nbytes) {
80 if( nbytes > cap ) {
81 if( ptr != NULL ) ascfree(ptr);
82 ptr = ASC_NEW_ARRAY(char,nbytes);
83 cap = nbytes;
84 }
85 } else {
86 if (ptr) ascfree(ptr);
87 ptr=NULL;
88 cap=0;
89 }
90 if ( cap >0) {
91 return(ptr);
92 } else {
93 return NULL;
94 }
95 }
96
97 #define rel_tmpalloc_array(nelts,type) \
98 ((nelts) > 0 ? (type *)tmpalloc((nelts)*sizeof(type)) : NULL)
99 /**
100 *** Creates an array of "nelts" objects, each with type "type".
101 **/
102
103 void relman_free_reused_mem(void)
104 {
105 /* rel_tmpalloc(0); */
106 RelationFindRoots(NULL,0,0,0,0,NULL,NULL,NULL);
107 }
108
109
110 #if REIMPLEMENT
111 boolean relman_is_linear( struct rel_relation *rel, var_filter_t *filter)
112 {
113 return (
114 exprman_is_linear(rel,rel_lhs(rel),filter) &&
115 exprman_is_linear(rel,rel_rhs(rel),filter)
116 );
117 }
118
119 real64 relman_linear_coef(struct rel_relation *rel, struct var_variable *var,
120 var_filter_t *filter)
121 {
122 return(
123 exprman_linear_coef(rel,rel_lhs(rel),var,filter) -
124 exprman_linear_coef(rel,rel_rhs(rel),var,filter)
125 );
126 }
127 #endif
128
129 #if KILL
130 void relman_decide_incidence( struct rel_relation *rel)
131 {
132 struct var_variable **list;
133 int c;
134
135 list = rel_incidence_list(rel);
136 for( c=0; list[c] != NULL; c++ )
137 var_set_incident(list[c],TRUE);
138 }
139 #endif
140
141 void relman_get_incidence(struct rel_relation *rel, var_filter_t *filter,
142 mtx_matrix_t mtx)
143 {
144 const struct var_variable **list;
145 mtx_coord_t nz;
146 int c,len;
147
148 assert(rel!=NULL && filter !=NULL && mtx != NULL);
149 nz.row = rel_sindex(rel);
150 len = rel_n_incidences(rel);
151
152 list = rel_incidence_list(rel);
153 for (c=0; c < len; c++) {
154 if( var_apply_filter(list[c],filter) ) {
155 nz.col = var_sindex(list[c]);
156 mtx_fill_org_value(mtx,&nz,1.0);
157 }
158 }
159 }
160
161 #ifdef RELOCATE_GB_NEEDED
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 = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n));
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 /** @TODO FIXME
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.
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 = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n)); /* 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 = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n));
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 #endif /* RELOCATE_GB_NEEDED */
372
373 real64 relman_eval(struct rel_relation *rel, int32 *status, int safe)
374 {
375 real64 res;
376 assert(status!=NULL && rel!=NULL);
377 if ( rel->type == e_rel_token ) {
378 if (!RelationCalcResidualBinary(
379 GetInstanceRelationOnly(IPTR(rel->instance)),&res)) {
380 *status = 1; /* calc_ok */
381 rel_set_residual(rel,res);
382 return res;
383 }
384 /* else we don't care -- go on to the old handling which
385 * is reasonably correct, if slow.
386 */
387 }
388 if( safe ) {
389 *status = (int32)RelationCalcResidualSafe(rel_instance(rel),&res);
390 if(*status)CONSOLE_DEBUG("Relation evaluation returned error");
391 /* CONSOLE_DEBUG("residual = %g",res); */
392 safe_error_to_stderr( (enum safe_err *)status );
393 /* always set the relation residual when using safe functions */
394 rel_set_residual(rel,res);
395 } else {
396 *status = RelationCalcResidual(rel_instance(rel),&res);
397 if ( *status ) {
398 /* an error occured */
399 res = 1.0e8;
400 } else {
401 rel_set_residual(rel,res);
402 }
403 }
404 /* flip the status flag: all values other than safe_ok become 0 */
405 *status = !(*status);
406 /* CONSOLE_DEBUG("returning %g",res); */
407 return res;
408
409 #if REIMPLEMENT /* all this is to be done on the compiler side, without
410 all the indirection idiocy. */
411 it may take some changes in the CalcResidual header to do so.
412 switch (rel->type) {
413 case e_token:
414 res = exprman_eval(rel,rel_lhs(rel)) - exprman_eval(rel,rel_rhs(rel));
415 break;
416 case e_opcode:
417 FPRINTF(stderr,"opcode relation processing not yet supported\n");
418 res = 1.0e08;
419 break;
420 case e_glassbox:
421 res = relman_glassbox_eval(rel);
422 break;
423 case e_blackbox:
424 res = ExtRel_Evaluate_LHS(rel) - ExtRel_Evaluate_RHS(rel);
425 break;
426 default:
427 FPRINTF(stderr,"unknown relation type in (relman_eval)\n");
428 res = 1.0e08;
429 break;
430 }
431 #endif
432
433 }
434 int32 relman_obj_direction(struct rel_relation *rel)
435 {
436 assert(rel!=NULL);
437 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
438 case e_minimize:
439 return -1;
440 case e_maximize:
441 return 1;
442 default:
443 return 0;
444 }
445 }
446
447 real64 relman_scale(struct rel_relation *rel)
448 {
449 real64 relnom;
450 assert(rel!=NULL);
451 relnom = CalcRelationNominal(rel_instance(rel));
452 if(relnom < 0.0000001) {
453 /* take care of small relnoms and relnom = 0 error returns */
454 relnom = 1.0;
455 }
456 rel_set_nominal(rel,relnom);
457 return relnom;
458 }
459
460 #if REIMPLEMENT /* compiler */
461 real64 relman_diff(struct rel_relation *rel, struct var_variable *var,
462 int safe)
463 {
464 /* FIX FIX FIX meaning kirk couldn't be botghered... */
465 real64 res = 0.0;
466 switch(rel->type) {
467 case e_glassbox:
468 case e_blackbox:
469 FPRINTF(stderr,"relman_diff not yet supported for black and glassbox\n");
470 return 1.0e08;
471 }
472 res -= exprman_diff(rel,rel_rhs(rel),var);
473 res += exprman_diff(rel,rel_lhs(rel),var);
474 return( res );
475 }
476 #endif
477
478 /* return 0 on success */
479 int relman_diff2(struct rel_relation *rel, var_filter_t *filter,
480 real64 *derivatives, int32 *variables,
481 int32 *count, int32 safe)
482 {
483 const struct var_variable **vlist=NULL;
484 real64 *gradient;
485 int32 len,c;
486 int status;
487
488 assert(rel!=NULL && filter!=NULL);
489 len = rel_n_incidences(rel);
490 vlist = rel_incidence_list(rel);
491
492 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
493 assert(gradient !=NULL);
494 *count = 0;
495 if(safe){
496 status =(int32)RelationCalcGradientSafe(rel_instance(rel),gradient);
497 safe_error_to_stderr( (enum safe_err *)&status );
498 /* always map when using safe functions */
499 for (c=0; c < len; c++) {
500 if (var_apply_filter(vlist[c],filter)) {
501 variables[*count] = var_sindex(vlist[c]);
502 derivatives[*count] = gradient[c];
503 /* CONSOLE_DEBUG("Var %d = %f",var_sindex(vlist[c]),gradient[c]); */
504 (*count)++;
505 }
506 }
507 /* CONSOLE_DEBUG("RETURNING (SAFE) calc_ok=%d",status); */
508 return status;
509 }else{
510 if((status=RelationCalcGradient(rel_instance(rel),gradient)) == 0) {
511 /* successful */
512 for (c=0; c < len; c++) {
513 if (var_apply_filter(vlist[c],filter)) {
514 variables[*count] = var_sindex(vlist[c]);
515 derivatives[*count] = gradient[c];
516 (*count)++;
517 }
518 }
519 }
520 /* SOLE_DEBUG("RETURNING (NON-SAFE) calc_ok=%d",status); */
521 return status;
522 }
523 }
524
525 int relman_diff_grad(struct rel_relation *rel, var_filter_t *filter,
526 real64 *derivatives, int32 *variables_master,
527 int32 *variables_solver, int32 *count, real64 *resid,
528 int32 safe)
529 {
530 const struct var_variable **vlist=NULL;
531 real64 *gradient;
532 int32 len,c;
533 int status;
534
535 assert(rel!=NULL && filter!=NULL);
536 len = rel_n_incidences(rel);
537 vlist = rel_incidence_list(rel);
538
539 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
540 assert(gradient !=NULL);
541 *count = 0;
542 if( safe ) {
543 /* CONSOLE_DEBUG("..."); */
544 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),
545 resid,gradient);
546 safe_error_to_stderr( (enum safe_err *)&status );
547 /* always map when using safe functions */
548 for (c=0; c < len; c++) {
549 if (var_apply_filter(vlist[c],filter)) {
550 variables_master[*count] = var_mindex(vlist[c]);
551 variables_solver[*count] = var_sindex(vlist[c]);
552 derivatives[*count] = gradient[c];
553 (*count)++;
554 }
555 }
556 }
557 else {
558 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient))== 0) {
559 /* successful */
560 for (c=0; c < len; c++) {
561 if (var_apply_filter(vlist[c],filter)) {
562 variables_master[*count] = var_mindex(vlist[c]);
563 variables_solver[*count] = var_sindex(vlist[c]);
564 derivatives[*count] = gradient[c];
565 (*count)++;
566 }
567 }
568 }
569 }
570
571 return !status; /* flip the status flag */
572 }
573
574 int32 relman_diff_harwell(struct rel_relation **rlist,
575 var_filter_t *vfilter, rel_filter_t *rfilter,
576 int32 rlen, int32 bias, int32 mORs,
577 real64 *avec, int32 *ivec, int32 *jvec)
578 {
579 const struct var_variable **vlist = NULL;
580 struct rel_relation *rel;
581 real64 residual, *resid;
582 real64 *gradient;
583 mtx_coord_t coord;
584 int32 len,c,r,k;
585 int32 errcnt;
586 enum safe_err status;
587
588 if(rlist == NULL || vfilter == NULL || rfilter == NULL || avec == NULL
589 || rlen < 0 || mORs >3 || mORs < 0 || bias <0 || bias > 1
590 ){
591 return 1;
592 }
593 resid = &residual;
594 errcnt = k = 0;
595
596 if(ivec==NULL || jvec==NULL){
597 /*_skip stuffing ivec,jvec */
598 for (r=0; r < rlen; r++) {
599 rel = rlist[r];
600 len = rel_n_incidences(rel);
601 vlist = rel_incidence_list(rel);
602 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
603 if (gradient == NULL) {
604 return 1;
605 }
606 /* CONSOLE_DEBUG("..."); */
607 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
608 safe_error_to_stderr(&status);
609 if (status) {
610 errcnt--;
611 }
612 for (c=0; c < len; c++) {
613 if (var_apply_filter(vlist[c],vfilter)) {
614 avec[k] = gradient[c];
615 k++;
616 }
617 }
618 }
619 }else{
620 for (r=0; r < rlen; r++) {
621 rel = rlist[r];
622 len = rel_n_incidences(rel);
623 vlist = rel_incidence_list(rel);
624 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
625 if (gradient == NULL) {
626 return 1;
627 }
628 /* CONSOLE_DEBUG("..."); */
629 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
630 safe_error_to_stderr(&status);
631 if (status) {
632 errcnt--;
633 }
634 if (mORs & 2) {
635 coord.row = rel_sindex(rel);
636 } else {
637 coord.row = rel_mindex(rel);
638 }
639 for (c=0; c < len; c++) {
640 if (var_apply_filter(vlist[c],vfilter)) {
641 if (mORs & 1) {
642 coord.col = var_sindex(vlist[c]);
643 } else {
644 coord.col = var_mindex(vlist[c]);
645 }
646 avec[k] = gradient[c];
647 ivec[k] = coord.row;
648 jvec[k] = coord.col;
649 k++;
650 }
651 }
652 }
653 }
654 return errcnt;
655 }
656
657 int32 relman_jacobian_count(struct rel_relation **rlist, int32 rlen,
658 var_filter_t *vfilter,
659 rel_filter_t *rfilter, int32 *max)
660 {
661 int32 len, result=0, row, count, c;
662 const struct var_variable **list;
663 struct rel_relation *rel;
664 *max = 0;
665
666 for( row = 0; row < rlen; row++ ) {
667 rel = rlist[row];
668 if (rel_apply_filter(rel,rfilter)) {
669 len = rel_n_incidences(rel);
670 list = rel_incidence_list(rel);
671 count = 0;
672 for (c=0; c < len; c++) {
673 if( var_apply_filter(list[c],vfilter) ) {
674 count++;
675 }
676 }
677 result += count;
678 *max = MAX(*max,count);
679 }
680 }
681 return result;
682 }
683
684 int relman_diffs(struct rel_relation *rel, var_filter_t *filter,
685 mtx_matrix_t mtx, real64 *resid, int safe)
686 {
687 const struct var_variable **vlist=NULL;
688 real64 *gradient;
689 int32 len,c;
690 mtx_coord_t coord;
691 int status;
692
693 assert(rel!=NULL && filter!=NULL && mtx != NULL);
694 len = rel_n_incidences(rel);
695 vlist = rel_incidence_list(rel);
696 coord.row = rel_sindex(rel);
697 assert(coord.row>=0 && coord.row < mtx_order(mtx));
698
699 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
700 assert(gradient !=NULL);
701 if( safe ) {
702 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
703 safe_error_to_stderr( (enum safe_err *)&status );
704 /* always map when using safe functions */
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 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 /* flip the status flag */
726 return !status;
727 }
728
729 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
730 real64 relman_diffs_orig( struct rel_relation *rel, var_filter_t *filter,
731 mtx_matrix_t mtx)
732 {
733 real64 res = 0.0;
734 int32 row;
735 row = mtx_org_to_row(mtx,rel_sindex(rel));
736 if (rel_extnodeinfo(rel)) {
737 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
738 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
739 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
740 return res;
741 } else {
742 res -= exprman_diffs(rel,rel_rhs(rel),filter,row,mtx);
743 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
744 res += exprman_diffs(rel,rel_lhs(rel),filter,row,mtx);
745 return(res);
746 }
747 }
748 #endif
749
750 boolean relman_calc_satisfied( struct rel_relation *rel, real64 tolerance)
751 {
752 real64 res;
753 res = rel_residual(rel);
754 if (!asc_finite(res)) {
755 rel_set_satisfied(rel,FALSE);
756 return( rel_satisfied(rel) );
757 }
758 if( res < 0.0 ) {
759 if( rel_less(rel) ) {
760 rel_set_satisfied(rel,TRUE);
761 return( rel_satisfied(rel) );
762 }
763 if( rel_greater(rel) ) {
764 rel_set_satisfied(rel,FALSE);
765 return( rel_satisfied(rel) );
766 }
767 rel_set_satisfied(rel,(-res <= tolerance ));
768 return( rel_satisfied(rel) );
769 }
770 if( res > 0.0 ) {
771 if( rel_greater(rel) ) {
772 rel_set_satisfied(rel,TRUE);
773 return( rel_satisfied(rel) );
774 }
775 if( rel_less(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 rel_set_satisfied(rel,rel_equal(rel)); /* strict >0 or <0 not satisfied */
783 return( rel_satisfied(rel) );
784 }
785
786 boolean relman_calc_satisfied_scaled(struct rel_relation *rel, real64 tolerance)
787 {
788 real64 res;
789 res = rel_residual(rel);
790
791 if( res < 0.0 )
792 if( rel_less(rel) )
793 rel_set_satisfied(rel,TRUE);
794 else if( rel_greater(rel) )
795 rel_set_satisfied(rel,FALSE);
796 else
797 rel_set_satisfied(rel,(-res <= tolerance * rel_nominal(rel)));
798 else if( res > 0.0 )
799 if( rel_greater(rel) )
800 rel_set_satisfied(rel,TRUE);
801 else if( rel_less(rel) )
802 rel_set_satisfied(rel,FALSE);
803 else
804 rel_set_satisfied(rel,(res <= tolerance * rel_nominal(rel)));
805 else
806 rel_set_satisfied(rel,rel_equal(rel));
807 return( rel_satisfied(rel) );
808 }
809
810 #if REIMPLEMENT
811 real64 *relman_directly_solve_new( struct rel_relation *rel,
812 struct var_variable *solvefor, int *able, int *nsolns,
813 real64 tolerance)
814 {
815 double *value;
816 if( rel_less(rel) || rel_greater(rel) || !rel_equal(rel) ||
817 rel_extnodeinfo(rel)) {
818 *able = FALSE;
819 *nsolns = 0;
820 return(NULL);
821 }
822 else if (rel->type == e_glassbox) {
823 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
824 return value;
825 }
826 else{
827 value =
828 exprman_directly_solve(rel,rel_lhs(rel),rel_rhs(rel),
829 solvefor,able,nsolns);
830 return value;
831 }
832 }
833
834 #else /* temporary */
835 real64 *relman_directly_solve_new( struct rel_relation *rel,
836 struct var_variable *solvefor, int *able, int *nsolns, real64 tolerance)
837 {
838 double *value;
839 if( rel_less(rel) ||
840 rel_greater(rel) ||
841 !rel_equal(rel) ||
842 rel_extnodeinfo(rel)) {
843 *able = FALSE;
844 *nsolns = 0;
845 return(NULL);
846 }
847 switch (rel->type ) {
848 #if 0 /* default */
849 case e_rel_glassbox:
850 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
851 return value;
852 #endif
853 case e_rel_token:
854 {
855 int nvars,n;
856 const struct var_variable **vlist;
857 unsigned long vindex; /* index to the compiler */
858 nvars = rel_n_incidences(rel);
859 vlist = rel_incidence_list(rel);
860 vindex = 0;
861 for (n=0;n <nvars; n++) {
862 if (vlist[n]==solvefor) {
863 vindex = n+1; /*compiler counts from 1 */
864 break;
865 }
866 }
867 value = RelationFindRoots(IPTR(rel_instance(rel)),
868 var_lower_bound(solvefor),
869 var_upper_bound(solvefor),
870 var_nominal(solvefor),
871 tolerance,
872 &(vindex),
873 able,nsolns);
874 return value;
875 }
876 default: /* e_rel_blackbox, glassbox */
877 *able = FALSE;
878 *nsolns = 0;
879 return(NULL);
880 }
881 }
882 #endif
883
884 char *dummyrelstring(slv_system_t sys, struct rel_relation *rel, int style)
885 {
886 char *result;
887 UNUSED_PARAMETER(sys);
888 UNUSED_PARAMETER(rel);
889 UNUSED_PARAMETER(style);
890 result = ASC_NEW_ARRAY(char,80);
891 sprintf(result,"relman_make_*string_*fix not implemented yet");
892 return result;
893 }
894
895 /*
896 * converst the relations vlist into an array of int indices suitable
897 * for relation write string (compiler side indexing from 1, remember).
898 * if mORs == TRUE master indices are used, else solver.
899 */
900 static
901 void relman_cookup_indices(struct RXNameData *rd,
902 struct rel_relation *rel,int mORs)
903 {
904 int nvars,n;
905 const struct var_variable **vlist;
906
907 nvars = rel_n_incidences(rel);
908 rd->indices = rel_tmpalloc_array(1+nvars,int);
909 if (rd->indices == NULL) {
910 return;
911 }
912 vlist = rel_incidence_list(rel);
913 if (mORs) {
914 for (n = 0; n < nvars; n++) {
915 rd->indices[n+1] = var_mindex(vlist[n]);
916 }
917 } else {
918 for (n = 0; n < nvars; n++) {
919 rd->indices[n+1] = var_sindex(vlist[n]);
920 }
921 }
922 }
923
924 char *relman_make_vstring_infix(slv_system_t sys,
925 struct rel_relation *rel, int style)
926 {
927 char *sbeg;
928 struct RXNameData rd = {"x",NULL,""};
929
930 if (style) {
931 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
932 NULL,NULL,relio_ascend,NULL);
933 } else {
934 /* need to cook up output indices */
935 relman_cookup_indices(&rd,rel,1);
936 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
937 (WRSNameFunc)RelationVarXName,
938 &rd,relio_ascend,NULL);
939 }
940 return(sbeg);
941 }
942
943 char *relman_make_vstring_postfix(slv_system_t sys,
944 struct rel_relation *rel, int style)
945 {
946 char *sbeg;
947
948 if (style) {
949 sbeg = WriteRelationPostfixString(rel_instance(rel),slv_instance(sys));
950 } else {
951 #if REIMPLEMENT
952 left = exprman_make_xstring_postfix(rel,sys,rel_lhs(rel));
953 right = exprman_make_xstring_postfix(rel,sys,rel_rhs(rel));
954 #else
955 sbeg = ASC_NEW_ARRAY(char,60);
956 if (sbeg==NULL) return sbeg;
957 sprintf(sbeg,"relman_make_xstring_postfix not reimplemented.");
958 #endif
959 }
960 return(sbeg);
961 }

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