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

Contents of /trunk/ascend/system/relman.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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