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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2351 - (show annotations) (download) (as text)
Thu Jan 6 02:02:41 2011 UTC (11 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 31104 byte(s)
Resolved memory leak with test/test solver_ipopt.formula.
Still some issues to overcome with regard to inactive relations, probably.
Some issue still occurring when attempting to run whole solver_ipopt suite in one go.
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 CONSOLE_DEBUG("Var %d = %g",var_sindex(vlist[c]),gradient[c]);
498 (*count)++;
499 }
500 }
501 /* CONSOLE_DEBUG("RETURNING (SAFE) calc_ok=%d",status); */
502 return status;
503 }else{
504 //CONSOLE_DEBUG("Derivative Type: Not SAFE");
505 if((status=RelationCalcGradient(rel_instance(rel),gradient)) == 0) {
506 /* successful */
507 for (c=0; c < len; c++) {
508 if (var_apply_filter(vlist[c],filter)) {
509 variables[*count] = var_sindex(vlist[c]);
510 derivatives[*count] = gradient[c];
511 (*count)++;
512 }
513 }
514 }
515 /* SOLE_DEBUG("RETURNING (NON-SAFE) calc_ok=%d",status); */
516 return status;
517 }
518 }
519
520
521 /* return 0 on success (derivatives, variables and count are output vars too) */
522 int relman_diff2_rev(struct rel_relation *rel, const var_filter_t *filter
523 ,real64 *derivatives, int32 *variables
524 ,int32 *count, int32 safe)
525 {
526 const struct var_variable **vlist=NULL;
527 real64 *gradient;
528 int32 len,c;
529 int status;
530 assert(rel!=NULL && filter!=NULL);
531 len = rel_n_incidences(rel);
532 // CONSOLE_DEBUG("In Function relman_diff2_rev");
533 vlist = rel_incidence_list(rel);
534
535 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
536 assert(gradient !=NULL);
537 *count = 0;
538 if(safe){
539 //CONSOLE_DEBUG("Derivative Type: Safe");
540 //PrintGradients(rel_instance(rel));
541 status =(int32)RelationCalcGradientRevSafe(rel_instance(rel),gradient);
542 safe_error_to_stderr( (enum safe_err *)&status );
543 /* always map when using safe functions */
544 for (c=0; c < len; c++) {
545 // CONSOLE_DEBUG("c = %d",c);
546 if (var_apply_filter(vlist[c],filter)) {
547 variables[*count] = var_sindex(vlist[c]);
548 derivatives[*count] = gradient[c];
549 // CONSOLE_DEBUG("Var %d = %g",var_sindex(vlist[c]),gradient[c]);
550 (*count)++;
551 }
552 }
553 // CONSOLE_DEBUG("RETURNING (SAFE) calc_ok=%d",status);
554 return status;
555 }else{
556 //CONSOLE_DEBUG("Derivative Type: Not SAFE");
557 if(
558 (status = (int32)RelationCalcGradientRev(rel_instance(rel),gradient))
559 == 0
560 ){
561 /* successful */
562 for (c=0; c < len; c++) {
563 if (var_apply_filter(vlist[c],filter)) {
564 variables[*count] = var_sindex(vlist[c]);
565 derivatives[*count] = gradient[c];
566 (*count)++;
567 }
568 }
569 }
570 /* CONSOLE_DEBUG("RETURNING (NON-SAFE) calc_ok=%d",status); */
571 return status;
572 }
573 }
574
575
576 /** ---------------------Hessian Calculations------------------------------ */
577 /* return 0 on success (derivatives, variables and count are output vars too) */
578 int relman_hess(struct rel_relation *rel, const var_filter_t *filter
579 ,hessian_mtx *hess_matrix,int32 *count,unsigned long max_dimension, int32 safe)
580 {
581 const struct var_variable **vlist=NULL;
582 hessian_mtx *matrix;
583 int32 len,i,j;
584 int status;
585
586 assert(rel!=NULL && filter!=NULL);
587
588 len = rel_n_incidences(rel);
589 #if 0
590 /* we can't apply this assertion because there may be unfiltered incidences making len
591 greater than max_dimension! */
592 asc_assert(len<=max_dimension); //Checking if Index is out of bounds
593 #endif
594
595 vlist = rel_incidence_list(rel);
596
597 // CONSOLE_DEBUG("IN FUNCTION relman_hess");
598
599 matrix = Hessian_Mtx_create(hess_matrix->access_type,len); // As Hessians may be (rarely) unsymmetrical
600 // type of Hessian matrix should be decided from
601 // type of relation
602 asc_assert(matrix !=NULL);
603 *count = 0;
604
605
606
607 if(safe){
608 status =(int32)RelationCalcHessianMtxSafe(rel_instance(rel),matrix,len);
609 safe_error_to_stderr( (enum safe_err *)&status );
610 /* always map when using safe functions */
611 for(i=0;i<len;i++){
612 if(var_apply_filter(vlist[i],filter)){
613 for(j=0;j<=i;j++){
614 if (var_apply_filter(vlist[j],filter)) {
615 Hessian_Mtx_set_element(hess_matrix,i,j,Hessian_Mtx_get_element(matrix,i,j));
616 (*count)++;
617 }
618 }
619 }
620 }
621 // CONSOLE_DEBUG("RETURNING (SAFE) calc_ok=%d",status);
622 }else{
623 if((status =(int32)RelationCalcHessianMtx(rel_instance(rel),matrix,len)) == 0) {
624
625 /* successful */
626 for(i=0;i<len;i++){
627 if(var_apply_filter(vlist[i],filter)){
628 for(j=0;j<=i;j++){
629 if (var_apply_filter(vlist[j],filter)) {
630 Hessian_Mtx_set_element(hess_matrix,i,j,Hessian_Mtx_get_element(matrix,i,j));
631 (*count)++;
632 }
633 }
634 }
635 }
636 }
637 }
638
639 Hessian_Mtx_destroy(matrix);
640
641 return status;
642 }
643
644 /* return 0 on success */
645 int relman_diff3(struct rel_relation *rel
646 , const var_filter_t *filter
647 , real64 *derivatives, struct var_variable **variables
648 , int32 *count, int32 safe
649 ){
650 struct var_variable **vlist=NULL;
651 real64 *gradient;
652 int32 len,c;
653 int status;
654
655 assert(rel!=NULL && filter!=NULL);
656 len = rel_n_incidences(rel);
657 vlist = (struct var_variable**)rel_incidence_list(rel);
658
659 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
660 assert(gradient !=NULL);
661 *count = 0;
662 if(safe){
663 #ifdef DIFF_DEBUG
664 CONSOLE_DEBUG("SAFE EVALUATION");
665 #endif
666 status =(int32)RelationCalcGradientSafe(rel_instance(rel),gradient);
667 safe_error_to_stderr( (enum safe_err *)&status );
668 /* always map when using safe functions */
669 for (c=0; c < len; c++) {
670 if (var_apply_filter(vlist[c],filter)) {
671 variables[*count] = vlist[c];
672 derivatives[*count] = gradient[c];
673 /* CONSOLE_DEBUG("Var %d = %f",var_sindex(vlist[c]),gradient[c]); */
674 (*count)++;
675 }
676 }
677 /* CONSOLE_DEBUG("RETURNING (SAFE) calc_ok=%d",status); */
678 return status;
679 }else{
680 #ifdef DIFF_DEBUG
681 CONSOLE_DEBUG("UNSAFE EVALUATION");
682 #endif
683 if((status=RelationCalcGradient(rel_instance(rel),gradient)) == 0) {
684 /* successful */
685 for (c=0; c < len; c++) {
686 if (var_apply_filter(vlist[c],filter)) {
687 variables[*count] = vlist[c];
688 derivatives[*count] = gradient[c];
689 (*count)++;
690 }
691 }
692 }
693 /* SOLE_DEBUG("RETURNING (NON-SAFE) calc_ok=%d",status); */
694 return status;
695 }
696 }
697
698
699 int relman_diff_grad(struct rel_relation *rel
700 , const var_filter_t *filter
701 , real64 *derivatives, int32 *variables_master
702 , int32 *variables_solver, int32 *count, real64 *resid
703 , int32 safe
704 ){
705 const struct var_variable **vlist=NULL;
706 real64 *gradient;
707 int32 len,c;
708 int status;
709
710 assert(rel!=NULL && filter!=NULL);
711 len = rel_n_incidences(rel);
712 vlist = rel_incidence_list(rel);
713
714 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
715 assert(gradient !=NULL);
716 *count = 0;
717 if( safe ) {
718 /* CONSOLE_DEBUG("..."); */
719 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),
720 resid,gradient);
721 safe_error_to_stderr( (enum safe_err *)&status );
722 /* always map when using safe functions */
723 for (c=0; c < len; c++) {
724 if (var_apply_filter(vlist[c],filter)) {
725 variables_master[*count] = var_mindex(vlist[c]);
726 variables_solver[*count] = var_sindex(vlist[c]);
727 derivatives[*count] = gradient[c];
728 (*count)++;
729 }
730 }
731 }
732 else {
733 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient))== 0) {
734 /* successful */
735 for (c=0; c < len; c++) {
736 if (var_apply_filter(vlist[c],filter)) {
737 variables_master[*count] = var_mindex(vlist[c]);
738 variables_solver[*count] = var_sindex(vlist[c]);
739 derivatives[*count] = gradient[c];
740 (*count)++;
741 }
742 }
743 }
744 }
745
746 return !status; /* flip the status flag */
747 }
748
749
750 int32 relman_diff_harwell(struct rel_relation **rlist
751 , var_filter_t *vfilter, rel_filter_t *rfilter
752 , int32 rlen, int32 bias, int32 mORs
753 , real64 *avec, int32 *ivec, int32 *jvec
754 ){
755 const struct var_variable **vlist = NULL;
756 struct rel_relation *rel;
757 real64 residual, *resid;
758 real64 *gradient;
759 mtx_coord_t coord;
760 int32 len,c,r,k;
761 int32 errcnt;
762 enum safe_err status;
763
764 if(rlist == NULL || vfilter == NULL || rfilter == NULL || avec == NULL
765 || rlen < 0 || mORs >3 || mORs < 0 || bias <0 || bias > 1
766 ){
767 return 1;
768 }
769 resid = &residual;
770 errcnt = k = 0;
771
772 if(ivec==NULL || jvec==NULL){
773 /*_skip stuffing ivec,jvec */
774 for (r=0; r < rlen; r++) {
775 rel = rlist[r];
776 len = rel_n_incidences(rel);
777 vlist = rel_incidence_list(rel);
778 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
779 if (gradient == NULL) {
780 return 1;
781 }
782 /* CONSOLE_DEBUG("..."); */
783 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
784 safe_error_to_stderr(&status);
785 if (status) {
786 errcnt--;
787 }
788 for (c=0; c < len; c++) {
789 if (var_apply_filter(vlist[c],vfilter)) {
790 avec[k] = gradient[c];
791 k++;
792 }
793 }
794 }
795 }else{
796 for (r=0; r < rlen; r++) {
797 rel = rlist[r];
798 len = rel_n_incidences(rel);
799 vlist = rel_incidence_list(rel);
800 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
801 if (gradient == NULL) {
802 return 1;
803 }
804 /* CONSOLE_DEBUG("..."); */
805 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
806 safe_error_to_stderr(&status);
807 if (status) {
808 errcnt--;
809 }
810 if (mORs & 2) {
811 coord.row = rel_sindex(rel);
812 }else{
813 coord.row = rel_mindex(rel);
814 }
815 for (c=0; c < len; c++) {
816 if (var_apply_filter(vlist[c],vfilter)) {
817 if (mORs & 1) {
818 coord.col = var_sindex(vlist[c]);
819 }else{
820 coord.col = var_mindex(vlist[c]);
821 }
822 avec[k] = gradient[c];
823 ivec[k] = coord.row;
824 jvec[k] = coord.col;
825 k++;
826 }
827 }
828 }
829 }
830 return errcnt;
831 }
832
833 int32 relman_jacobian_count(struct rel_relation **rlist, int32 rlen
834 , var_filter_t *vfilter
835 , rel_filter_t *rfilter, int32 *max
836 ){
837 int32 len, result=0, row, count, c;
838 const struct var_variable **list;
839 struct rel_relation *rel;
840 *max = 0;
841
842 for( row = 0; row < rlen; row++ ) {
843 rel = rlist[row];
844 if (rel_apply_filter(rel,rfilter)) {
845 len = rel_n_incidences(rel);
846 list = rel_incidence_list(rel);
847 count = 0;
848 for (c=0; c < len; c++) {
849 if( var_apply_filter(list[c],vfilter) ) {
850 count++;
851 }
852 }
853
854 result += count;
855 *max = MAX(*max,count);
856 }
857 }
858 return result;
859 }
860
861
862 static int AllVariables(struct Instance *i){
863 return TRUE;
864 }
865
866 int32 relman_hessian_count(struct rel_relation **rlist, int32 rlen
867 , var_filter_t *vfilter
868 , rel_filter_t *rfilter, int32 *max
869 ){
870
871 struct relation *r1, *r2;
872
873 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"not implemented...");
874
875 r1 = rel_instance(rlist[0]);
876
877 CONSOLE_DEBUG("r1 = %p",r1);
878 r2 = RelDerivative(r1,0,&AllVariables);
879
880 return 0;
881 }
882
883
884 int relman_diffs(struct rel_relation *rel
885 , const var_filter_t *filter
886 , mtx_matrix_t mtx, real64 *resid, int safe
887 ){
888 const struct var_variable **vlist=NULL;
889 real64 *gradient;
890 int32 len,c;
891 mtx_coord_t coord;
892 int status;
893
894 assert(rel!=NULL && filter!=NULL && mtx != NULL);
895 len = rel_n_incidences(rel);
896 vlist = rel_incidence_list(rel);
897 coord.row = rel_sindex(rel);
898 assert(coord.row>=0 && coord.row < mtx_order(mtx));
899
900 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
901 assert(gradient !=NULL);
902 if( safe ) {
903 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
904 safe_error_to_stderr( (enum safe_err *)&status );
905 /* always map when using safe functions */
906 for (c=0; c < len; c++) {
907 if (var_apply_filter(vlist[c],filter)) {
908 coord.col = var_sindex(vlist[c]);
909 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
910 mtx_fill_org_value(mtx,&coord,gradient[c]);
911 }
912 }
913 }
914 else {
915 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient)) == 0) {
916 /* successful */
917 for (c=0; c < len; c++) {
918 if (var_apply_filter(vlist[c],filter)) {
919 coord.col = var_sindex(vlist[c]);
920 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
921 mtx_fill_org_value(mtx,&coord,gradient[c]);
922 }
923 }
924 }
925 }
926 /* flip the status flag */
927 return !status;
928 }
929
930
931 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
932 real64 relman_diffs_orig( struct rel_relation *rel, var_filter_t *filter
933 ,mtx_matrix_t mtx
934 ){
935 real64 res = 0.0;
936 int32 row;
937 row = mtx_org_to_row(mtx,rel_sindex(rel));
938 if (rel_extnodeinfo(rel)) {
939 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
940 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
941 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
942 return res;
943 }else{
944 res -= exprman_diffs(rel,rel_rhs(rel),filter,row,mtx);
945 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
946 res += exprman_diffs(rel,rel_lhs(rel),filter,row,mtx);
947 return(res);
948 }
949 }
950 #endif
951
952
953 boolean relman_calc_satisfied( struct rel_relation *rel, real64 tolerance){
954 real64 res;
955 res = rel_residual(rel);
956 if (!asc_finite(res)) {
957 rel_set_satisfied(rel,FALSE);
958 return( rel_satisfied(rel) );
959 }
960 if( res < 0.0 ) {
961 if( rel_less(rel) ) {
962 rel_set_satisfied(rel,TRUE);
963 return( rel_satisfied(rel) );
964 }
965 if( rel_greater(rel) ) {
966 rel_set_satisfied(rel,FALSE);
967 return( rel_satisfied(rel) );
968 }
969 rel_set_satisfied(rel,(-res <= tolerance ));
970 return( rel_satisfied(rel) );
971 }
972 if( res > 0.0 ) {
973 if( rel_greater(rel) ) {
974 rel_set_satisfied(rel,TRUE);
975 return( rel_satisfied(rel) );
976 }
977 if( rel_less(rel) ) {
978 rel_set_satisfied(rel,FALSE);
979 return( rel_satisfied(rel) );
980 }
981 rel_set_satisfied(rel,(res <= tolerance ));
982 return( rel_satisfied(rel) );
983 }
984 rel_set_satisfied(rel,rel_equal(rel)); /* strict >0 or <0 not satisfied */
985 return( rel_satisfied(rel) );
986 }
987
988 boolean relman_calc_satisfied_scaled(struct rel_relation *rel, real64 tolerance)
989 {
990 real64 res;
991 res = rel_residual(rel);
992
993 if( res < 0.0 )
994 if( rel_less(rel) )
995 rel_set_satisfied(rel,TRUE);
996 else if( rel_greater(rel) )
997 rel_set_satisfied(rel,FALSE);
998 else
999 rel_set_satisfied(rel,(-res <= tolerance * rel_nominal(rel)));
1000 else if( res > 0.0 )
1001 if( rel_greater(rel) )
1002 rel_set_satisfied(rel,TRUE);
1003 else if( rel_less(rel) )
1004 rel_set_satisfied(rel,FALSE);
1005 else
1006 rel_set_satisfied(rel,(res <= tolerance * rel_nominal(rel)));
1007 else
1008 rel_set_satisfied(rel,rel_equal(rel));
1009 return( rel_satisfied(rel) );
1010 }
1011
1012 /* temporary */
1013 real64 *relman_directly_solve_new( struct rel_relation *rel,
1014 struct var_variable *solvefor, int *able, int *nsolns, real64 tolerance
1015 ){
1016 double *value;
1017 if(rel_less(rel) ||
1018 rel_greater(rel) ||
1019 !rel_equal(rel) ||
1020 rel_extnodeinfo(rel)
1021 ){
1022 /* fall through; not able to directly solve */
1023 }else{
1024 switch (rel->type ) {
1025 case e_rel_token:
1026 {
1027 int nvars,n;
1028 const struct var_variable **vlist;
1029 unsigned long vindex; /* index to the compiler */
1030 nvars = rel_n_incidences(rel);
1031 vlist = rel_incidence_list(rel);
1032 vindex = 0;
1033 for (n=0; n < nvars; n++) {
1034 if (vlist[n]==solvefor) {
1035 vindex = n+1; /*compiler counts from 1 */
1036 break;
1037 }
1038 }
1039 value = RelationFindRoots(IPTR(rel_instance(rel))
1040 , var_lower_bound(solvefor)
1041 , var_upper_bound(solvefor)
1042 , var_nominal(solvefor)
1043 , tolerance, &(vindex), able, nsolns
1044 );
1045 return value;
1046 }
1047 case e_rel_blackbox:
1048 #ifdef DSOLVE_DEBUG
1049 CONSOLE_DEBUG("Attempting direct solve of blackbox");
1050 #endif
1051 return blackbox_dsolve(
1052 IPTR(rel_instance(rel))
1053 ,IPTR(var_instance(solvefor))
1054 ,able,nsolns
1055 );
1056 /*
1057 we *could* directly solve if the solvefor variable is the output
1058 of this particular relation, but we don't support that at this stage
1059 */
1060 break;
1061 #if 0
1062 case e_rel_glassbox:
1063 /* I think glassbox functionality might be dead at the moment? */
1064 break;
1065 #endif
1066 default:
1067 /* anything else? */
1068 break;
1069 }
1070 }
1071 *able = FALSE;
1072 *nsolns = 0;
1073 return(NULL);
1074 }
1075
1076
1077 char *dummyrelstring(slv_system_t sys, struct rel_relation *rel, int style){
1078 char *result;
1079 UNUSED_PARAMETER(sys);
1080 UNUSED_PARAMETER(rel);
1081 UNUSED_PARAMETER(style);
1082 result = ASC_NEW_ARRAY(char,80);
1083 sprintf(result,"relman_make_*string_*fix not implemented yet");
1084 return result;
1085 }
1086
1087 /*
1088 * converst the relations vlist into an array of int indices suitable
1089 * for relation write string (compiler side indexing from 1, remember).
1090 * if mORs == TRUE master indices are used, else solver.
1091 */
1092 static
1093 void relman_cookup_indices(struct RXNameData *rd
1094 ,struct rel_relation *rel,int mORs
1095 ){
1096 int nvars,n;
1097 const struct var_variable **vlist;
1098
1099 nvars = rel_n_incidences(rel);
1100 rd->indices = rel_tmpalloc_array(1+nvars,int);
1101 if (rd->indices == NULL) {
1102 return;
1103 }
1104 vlist = rel_incidence_list(rel);
1105 if (mORs) {
1106 for (n = 0; n < nvars; n++) {
1107 rd->indices[n+1] = var_mindex(vlist[n]);
1108 }
1109 }else{
1110 for (n = 0; n < nvars; n++) {
1111 rd->indices[n+1] = var_sindex(vlist[n]);
1112 }
1113 }
1114 }
1115
1116
1117 char *relman_make_vstring_infix(slv_system_t sys
1118 ,struct rel_relation *rel, int style
1119 ){
1120 char *sbeg;
1121 struct RXNameData rd = {"x",NULL,""};
1122
1123 if (style) {
1124 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
1125 NULL,NULL,relio_ascend,NULL);
1126 }else{
1127 /* need to cook up output indices */
1128 relman_cookup_indices(&rd,rel,1);
1129 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
1130 (WRSNameFunc)RelationVarXName,
1131 &rd,relio_ascend,NULL);
1132 }
1133 return(sbeg);
1134 }
1135
1136
1137 char *relman_make_vstring_postfix(slv_system_t sys,
1138 struct rel_relation *rel, int style
1139 ){
1140 char *sbeg;
1141
1142 if (style) {
1143 sbeg = WriteRelationPostfixString(rel_instance(rel),slv_instance(sys));
1144 }else{
1145 #if REIMPLEMENT
1146 left = exprman_make_xstring_postfix(rel,sys,rel_lhs(rel));
1147 right = exprman_make_xstring_postfix(rel,sys,rel_rhs(rel));
1148 #else
1149 sbeg = ASC_NEW_ARRAY(char,60);
1150 if (sbeg==NULL) return sbeg;
1151 sprintf(sbeg,"relman_make_xstring_postfix not reimplemented.");
1152 #endif
1153 }
1154 return(sbeg);
1155 }

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