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

Contents of /trunk/base/generic/system/relman.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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