/[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 979 - (show annotations) (download) (as text)
Wed Dec 20 14:34:16 2006 UTC (13 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 26507 byte(s)
Added simplified ASC_PANIC call that uses var-args, added throughout relation_util.c.
Fixed var_filter_t stuff in djex and fvex.
More assertions in integrator.c
Added output of initial state from lsode.c (hoping that's a good idea?)
Fixed output code from relman_diff2.
Added asc_panic_nofunc for non var-arg CPPs.
Disabled -O3 flag in building C++ API
Added __getitem__ and __getattr__ methods in Simuluation for simplified python syntax (eg M.x instead M.sim.x)
Integrator::analyse throws exceptions on error now.

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

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