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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 909 - (show annotations) (download) (as text)
Thu Oct 26 12:44:41 2006 UTC (17 years, 9 months ago) by johnpye
File MIME type: text/x-csrc
File size: 26506 byte(s)
Added finite-difference evaluation of gradients in blackboxes.
Some work on J*v evaluation with IDA.
1 /*
2 * SLV: Ascend Nonlinear Solver
3 * by Karl Michael Westerberg
4 * Created: 2/6/90
5 * Version: $Revision: 1.44 $
6 * Version control file: $RCSfile: relman.c,v $
7 * Date last modified: $Date: 1998/04/23 23:56:22 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 *
16 * The SLV solver is free software; you can redistribute
17 * it and/or modify it under the terms of the GNU General Public License as
18 * published by the Free Software Foundation; either version 2 of the
19 * License, or (at your option) any later version.
20 *
21 * The SLV solver is distributed in hope that it will be
22 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24 * General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License along with
27 * the program; if not, write to the Free Software Foundation, Inc., 675
28 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
29 * COPYING is found in ../compiler.
30 */
31
32 #include <math.h>
33 #include <utilities/ascConfig.h>
34 #include <compiler/instance_enum.h>
35 #include <compiler/fractions.h>
36 #include <compiler/compiler.h>
37 #include <utilities/ascMalloc.h>
38 #include <compiler/functype.h>
39 #include <compiler/safe.h>
40 #include <general/list.h>
41 #include <compiler/extfunc.h>
42 #include <compiler/dimen.h>
43 #include <compiler/expr_types.h>
44 #include <compiler/find.h>
45 #include <compiler/atomvalue.h>
46 #include <compiler/mathinst.h>
47 #include <compiler/relation_type.h>
48 #include <compiler/rel_blackbox.h>
49 #include <compiler/vlist.h>
50 #include <compiler/relation.h>
51 #include <compiler/relation_util.h>
52 #include <compiler/relation_io.h>
53 #define _SLV_SERVER_C_SEEN_
54 #include "mtx.h"
55 #include "slv_types.h"
56 #include "var.h"
57 #include "rel.h"
58 #include "discrete.h"
59 #include "conditional.h"
60 #include "logrel.h"
61 #include "bnd.h"
62 #include "relman.h"
63 #include <compiler/rootfind.h>
64 #include "slv_server.h"
65 #include <general/mathmacros.h>
66
67 #define IPTR(i) ((struct Instance *)(i))
68
69 #define KILL 0 /* compile dead code if kill = 1 */
70 #define REIMPLEMENT 0 /* code that needs to be reimplemented */
71
72 static POINTER rel_tmpalloc( int nbytes)
73 /**
74 *** Temporarily allocates a given number of bytes. The memory need
75 *** not be freed, but the next call to this function will reuse the
76 *** previous allocation. Memory returned will NOT be zeroed.
77 *** Calling with nbytes==0 will free any memory allocated.
78 **/
79 {
80 static char *ptr = NULL;
81 static int cap = 0;
82
83 if (nbytes) {
84 if( nbytes > cap ) {
85 if( ptr != NULL ) ascfree(ptr);
86 ptr = ASC_NEW_ARRAY(char,nbytes);
87 cap = nbytes;
88 }
89 } else {
90 if (ptr) ascfree(ptr);
91 ptr=NULL;
92 cap=0;
93 }
94 if ( cap >0) {
95 return(ptr);
96 } else {
97 return NULL;
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 void relman_free_reused_mem(void)
108 {
109 /* rel_tmpalloc(0); */
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 {
117 return (
118 exprman_is_linear(rel,rel_lhs(rel),filter) &&
119 exprman_is_linear(rel,rel_rhs(rel),filter)
120 );
121 }
122
123 real64 relman_linear_coef(struct rel_relation *rel, struct var_variable *var,
124 var_filter_t *filter)
125 {
126 return(
127 exprman_linear_coef(rel,rel_lhs(rel),var,filter) -
128 exprman_linear_coef(rel,rel_rhs(rel),var,filter)
129 );
130 }
131 #endif
132
133 #if KILL
134 void relman_decide_incidence( struct rel_relation *rel)
135 {
136 struct var_variable **list;
137 int c;
138
139 list = rel_incidence_list(rel);
140 for( c=0; list[c] != NULL; c++ )
141 var_set_incident(list[c],TRUE);
142 }
143 #endif
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 #ifdef RELOCATE_GB_NEEDED
166 /*
167 *********************************************************************
168 * Code to deal with glassbox relations processing.
169 *********************************************************************
170 */
171
172 static double dsolve_scratch = 0.0; /* some workspace */
173 #define DSOLVE_TOLERANCE 1.0e-08 /* no longer needed */
174
175 static
176 real64 *relman_glassbox_dsolve(struct rel_relation *rel,
177 struct var_variable *solvefor,
178 int *able,
179 int *nsolns,
180 real64 tolerance)
181 {
182 int n,m,j,dummy = 0;
183 int mode, result;
184 int index,solve_for_index = -1;
185 struct Instance *var;
186 CONST struct relation *cmplr_reln;
187 enum Expr_enum type;
188 CONST struct gl_list_t *vlist;
189 double *f, *x, *g, value;
190 double lower, upper, nominal;
191 ExtEvalFunc **table;
192
193 cmplr_reln = GetInstanceRelation(rel_instance(rel),&type);
194 assert(type==e_glassbox);
195
196 vlist = RelationVarList(cmplr_reln);
197 m = rel_sindex(rel);
198 n = (int)gl_length(vlist);
199 f = ASC_NEW_ARRAY_CLEAR(double,(1 + 2*n));
200 x = &f[1];
201 g = &f[n+1];
202
203 /*
204 * Load x vector, and get the index into the x[vector]
205 * of the variable that we are solving for.
206 */
207 for (j=0;j<n;j++) {
208 var = (struct Instance *)gl_fetch(vlist,(unsigned long)(j+1));
209 x[j] = RealAtomValue(var);
210 if (var_instance(solvefor)== var)
211 solve_for_index = j;
212 }
213
214 /*
215 * Set up bounds and tolerance.
216 */
217 lower = var_lower_bound(solvefor);
218 upper = var_upper_bound(solvefor);
219 nominal = var_nominal(solvefor);
220 /*tolerance = DSOLVE_TOLERANCE; now passed in. */
221
222 /*
223 * Get the evaluation function.
224 */
225 table = GetValueJumpTable(GlassBoxExtFunc(cmplr_reln));
226 index = GlassBoxRelIndex(cmplr_reln);
227
228 /** @TODO FIXME
229 Call the rootfinder. A note of CAUTION:
230 zbrent is going to call the evaluation function.
231 It expects, that the results of this evaluation will be
232 written to f[m]. GlassBox relations however always
233 write to slot f[0]. To keep the world happy we send in
234 a dummy = 0. This needs to be made cleaner.
235 */
236 value = zbrent(table[index], &lower, &upper,
237 &mode, &dummy, &solve_for_index,
238 x, x , f, g, &tolerance, &result);
239 if (result==0) {
240 dsolve_scratch = value;
241 ascfree((char *)f);
242 *able = TRUE;
243 *nsolns = 1;
244 return &dsolve_scratch;
245 }
246 else{
247 dsolve_scratch = 0.0;
248 ascfree((char *)f);
249 *able = FALSE;
250 *nsolns = 0;
251 return NULL;
252 }
253 }
254
255
256 #ifdef THIS_IS_AN_UNUSED_FUNCTION
257 static
258 real64 relman_glassbox_eval(struct rel_relation *rel)
259 {
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 real64 relman_eval(struct rel_relation *rel, int32 *status, int safe)
378 {
379 real64 res;
380 assert(status!=NULL && rel!=NULL);
381 if ( rel->type == e_rel_token ) {
382 if (!RelationCalcResidualBinary(
383 GetInstanceRelationOnly(IPTR(rel->instance)),&res)) {
384 *status = 1; /* calc_ok */
385 rel_set_residual(rel,res);
386 return res;
387 }
388 /* else we don't care -- go on to the old handling which
389 * is reasonably correct, if slow.
390 */
391 }
392 if( safe ) {
393 *status = (int32)RelationCalcResidualSafe(rel_instance(rel),&res);
394 /* CONSOLE_DEBUG("residual = %g",res); */
395 safe_error_to_stderr( (enum safe_err *)status );
396 /* always set the relation residual when using safe functions */
397 rel_set_residual(rel,res);
398 } else {
399 *status = RelationCalcResidual(rel_instance(rel),&res);
400 if ( *status ) {
401 /* an error occured */
402 res = 1.0e8;
403 } else {
404 rel_set_residual(rel,res);
405 }
406 }
407 /* flip the status flag: all values other than safe_ok become 0 */
408 *status = !(*status);
409 /* CONSOLE_DEBUG("returning %g",res); */
410 return res;
411
412 #if REIMPLEMENT /* all this is to be done on the compiler side, without
413 all the indirection idiocy. */
414 it may take some changes in the CalcResidual header to do so.
415 switch (rel->type) {
416 case e_token:
417 res = exprman_eval(rel,rel_lhs(rel)) - exprman_eval(rel,rel_rhs(rel));
418 break;
419 case e_opcode:
420 FPRINTF(stderr,"opcode relation processing not yet supported\n");
421 res = 1.0e08;
422 break;
423 case e_glassbox:
424 res = relman_glassbox_eval(rel);
425 break;
426 case e_blackbox:
427 res = ExtRel_Evaluate_LHS(rel) - ExtRel_Evaluate_RHS(rel);
428 break;
429 default:
430 FPRINTF(stderr,"unknown relation type in (relman_eval)\n");
431 res = 1.0e08;
432 break;
433 }
434 #endif
435
436 }
437 int32 relman_obj_direction(struct rel_relation *rel)
438 {
439 assert(rel!=NULL);
440 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
441 case e_minimize:
442 return -1;
443 case e_maximize:
444 return 1;
445 default:
446 return 0;
447 }
448 }
449
450 real64 relman_scale(struct rel_relation *rel)
451 {
452 real64 relnom;
453 assert(rel!=NULL);
454 relnom = CalcRelationNominal(rel_instance(rel));
455 if(relnom < 0.0000001) {
456 /* take care of small relnoms and relnom = 0 error returns */
457 relnom = 1.0;
458 }
459 rel_set_nominal(rel,relnom);
460 return relnom;
461 }
462
463 #if REIMPLEMENT /* compiler */
464 real64 relman_diff(struct rel_relation *rel, struct var_variable *var,
465 int safe)
466 {
467 /* FIX FIX FIX meaning kirk couldn't be botghered... */
468 real64 res = 0.0;
469 switch(rel->type) {
470 case e_glassbox:
471 case e_blackbox:
472 FPRINTF(stderr,"relman_diff not yet supported for black and glassbox\n");
473 return 1.0e08;
474 }
475 res -= exprman_diff(rel,rel_rhs(rel),var);
476 res += exprman_diff(rel,rel_lhs(rel),var);
477 return( res );
478 }
479 #endif
480
481 int relman_diff2(struct rel_relation *rel, var_filter_t *filter,
482 real64 *derivatives, int32 *variables,
483 int32 *count, int32 safe)
484 {
485 const struct var_variable **vlist=NULL;
486 real64 *gradient;
487 int32 len,c;
488 int status;
489
490 assert(rel!=NULL && filter!=NULL);
491 len = rel_n_incidences(rel);
492 vlist = rel_incidence_list(rel);
493
494 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
495 assert(gradient !=NULL);
496 *count = 0;
497 if(safe){
498 status =(int32)RelationCalcGradientSafe(rel_instance(rel),gradient);
499 safe_error_to_stderr( (enum safe_err *)&status );
500 /* always map when using safe functions */
501 for (c=0; c < len; c++) {
502 if (var_apply_filter(vlist[c],filter)) {
503 variables[*count] = var_sindex(vlist[c]);
504 derivatives[*count] = gradient[c];
505 CONSOLE_DEBUG("Var %d = %f",var_sindex(vlist[c]),gradient[c]);
506 (*count)++;
507 }
508 }
509 return status;
510 }else{
511 if((status=RelationCalcGradient(rel_instance(rel),gradient)) == 0) {
512 /* successful */
513 for (c=0; c < len; c++) {
514 if (var_apply_filter(vlist[c],filter)) {
515 variables[*count] = var_sindex(vlist[c]);
516 derivatives[*count] = gradient[c];
517 (*count)++;
518 }
519 }
520 }
521 return !status;
522 }
523 }
524
525 int relman_diff_grad(struct rel_relation *rel, var_filter_t *filter,
526 real64 *derivatives, int32 *variables_master,
527 int32 *variables_solver, int32 *count, real64 *resid,
528 int32 safe)
529 {
530 const struct var_variable **vlist=NULL;
531 real64 *gradient;
532 int32 len,c;
533 int status;
534
535 assert(rel!=NULL && filter!=NULL);
536 len = rel_n_incidences(rel);
537 vlist = rel_incidence_list(rel);
538
539 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
540 assert(gradient !=NULL);
541 *count = 0;
542 if( safe ) {
543 /* CONSOLE_DEBUG("..."); */
544 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),
545 resid,gradient);
546 safe_error_to_stderr( (enum safe_err *)&status );
547 /* always map when using safe functions */
548 for (c=0; c < len; c++) {
549 if (var_apply_filter(vlist[c],filter)) {
550 variables_master[*count] = var_mindex(vlist[c]);
551 variables_solver[*count] = var_sindex(vlist[c]);
552 derivatives[*count] = gradient[c];
553 (*count)++;
554 }
555 }
556 }
557 else {
558 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient))== 0) {
559 /* successful */
560 for (c=0; c < len; c++) {
561 if (var_apply_filter(vlist[c],filter)) {
562 variables_master[*count] = var_mindex(vlist[c]);
563 variables_solver[*count] = var_sindex(vlist[c]);
564 derivatives[*count] = gradient[c];
565 (*count)++;
566 }
567 }
568 }
569 }
570
571 return !status; /* flip the status flag */
572 }
573
574 int32 relman_diff_harwell(struct rel_relation **rlist,
575 var_filter_t *vfilter, rel_filter_t *rfilter,
576 int32 rlen, int32 bias, int32 mORs,
577 real64 *avec, int32 *ivec, int32 *jvec)
578 {
579 const struct var_variable **vlist = NULL;
580 struct rel_relation *rel;
581 real64 residual, *resid;
582 real64 *gradient;
583 mtx_coord_t coord;
584 int32 len,c,r,k;
585 int32 errcnt;
586 enum safe_err status;
587
588 if(rlist == NULL || vfilter == NULL || rfilter == NULL || avec == NULL
589 || rlen < 0 || mORs >3 || mORs < 0 || bias <0 || bias > 1
590 ){
591 return 1;
592 }
593 resid = &residual;
594 errcnt = k = 0;
595
596 if(ivec==NULL || jvec==NULL){
597 /*_skip stuffing ivec,jvec */
598 for (r=0; r < rlen; r++) {
599 rel = rlist[r];
600 len = rel_n_incidences(rel);
601 vlist = rel_incidence_list(rel);
602 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
603 if (gradient == NULL) {
604 return 1;
605 }
606 /* CONSOLE_DEBUG("..."); */
607 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
608 safe_error_to_stderr(&status);
609 if (status) {
610 errcnt--;
611 }
612 for (c=0; c < len; c++) {
613 if (var_apply_filter(vlist[c],vfilter)) {
614 avec[k] = gradient[c];
615 k++;
616 }
617 }
618 }
619 }else{
620 for (r=0; r < rlen; r++) {
621 rel = rlist[r];
622 len = rel_n_incidences(rel);
623 vlist = rel_incidence_list(rel);
624 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
625 if (gradient == NULL) {
626 return 1;
627 }
628 /* CONSOLE_DEBUG("..."); */
629 status = RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
630 safe_error_to_stderr(&status);
631 if (status) {
632 errcnt--;
633 }
634 if (mORs & 2) {
635 coord.row = rel_sindex(rel);
636 } else {
637 coord.row = rel_mindex(rel);
638 }
639 for (c=0; c < len; c++) {
640 if (var_apply_filter(vlist[c],vfilter)) {
641 if (mORs & 1) {
642 coord.col = var_sindex(vlist[c]);
643 } else {
644 coord.col = var_mindex(vlist[c]);
645 }
646 avec[k] = gradient[c];
647 ivec[k] = coord.row;
648 jvec[k] = coord.col;
649 k++;
650 }
651 }
652 }
653 }
654 return errcnt;
655 }
656
657 int32 relman_jacobian_count(struct rel_relation **rlist, int32 rlen,
658 var_filter_t *vfilter,
659 rel_filter_t *rfilter, int32 *max)
660 {
661 int32 len, result=0, row, count, c;
662 const struct var_variable **list;
663 struct rel_relation *rel;
664 *max = 0;
665
666 for( row = 0; row < rlen; row++ ) {
667 rel = rlist[row];
668 if (rel_apply_filter(rel,rfilter)) {
669 len = rel_n_incidences(rel);
670 list = rel_incidence_list(rel);
671 count = 0;
672 for (c=0; c < len; c++) {
673 if( var_apply_filter(list[c],vfilter) ) {
674 count++;
675 }
676 }
677 result += count;
678 *max = MAX(*max,count);
679 }
680 }
681 return result;
682 }
683
684 int relman_diffs(struct rel_relation *rel, var_filter_t *filter,
685 mtx_matrix_t mtx, real64 *resid, int safe)
686 {
687 const struct var_variable **vlist=NULL;
688 real64 *gradient;
689 int32 len,c;
690 mtx_coord_t coord;
691 int status;
692
693 assert(rel!=NULL && filter!=NULL && mtx != NULL);
694 len = rel_n_incidences(rel);
695 vlist = rel_incidence_list(rel);
696 coord.row = rel_sindex(rel);
697 assert(coord.row>=0 && coord.row < mtx_order(mtx));
698
699 gradient = (real64 *)rel_tmpalloc(len*sizeof(real64));
700 assert(gradient !=NULL);
701 if( safe ) {
702 status =(int32)RelationCalcResidGradSafe(rel_instance(rel),resid,gradient);
703 safe_error_to_stderr( (enum safe_err *)&status );
704 /* always map when using safe functions */
705 for (c=0; c < len; c++) {
706 if (var_apply_filter(vlist[c],filter)) {
707 coord.col = var_sindex(vlist[c]);
708 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
709 mtx_fill_org_value(mtx,&coord,gradient[c]);
710 }
711 }
712 }
713 else {
714 if((status=RelationCalcResidGrad(rel_instance(rel),resid,gradient)) == 0) {
715 /* successful */
716 for (c=0; c < len; c++) {
717 if (var_apply_filter(vlist[c],filter)) {
718 coord.col = var_sindex(vlist[c]);
719 assert(coord.col >= 0 && coord.col < mtx_order(mtx));
720 mtx_fill_org_value(mtx,&coord,gradient[c]);
721 }
722 }
723 }
724 }
725 /* flip the status flag */
726 return !status;
727 }
728
729 #if REIMPLEMENT /* this needs to be reimplemented in the compiler */
730 real64 relman_diffs_orig( struct rel_relation *rel, var_filter_t *filter,
731 mtx_matrix_t mtx)
732 {
733 real64 res = 0.0;
734 int32 row;
735 row = mtx_org_to_row(mtx,rel_sindex(rel));
736 if (rel_extnodeinfo(rel)) {
737 res -= ExtRel_Diffs_RHS(rel,filter,row,mtx);
738 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
739 res += ExtRel_Diffs_LHS(rel,filter,row,mtx);
740 return res;
741 } else {
742 res -= exprman_diffs(rel,rel_rhs(rel),filter,row,mtx);
743 mtx_mult_row(mtx,row,-1.0,mtx_ALL_COLS);
744 res += exprman_diffs(rel,rel_lhs(rel),filter,row,mtx);
745 return(res);
746 }
747 }
748 #endif
749
750 boolean relman_calc_satisfied( struct rel_relation *rel, real64 tolerance)
751 {
752 real64 res;
753 res = rel_residual(rel);
754 if (!finite(res)) {
755 rel_set_satisfied(rel,FALSE);
756 return( rel_satisfied(rel) );
757 }
758 if( res < 0.0 ) {
759 if( rel_less(rel) ) {
760 rel_set_satisfied(rel,TRUE);
761 return( rel_satisfied(rel) );
762 }
763 if( rel_greater(rel) ) {
764 rel_set_satisfied(rel,FALSE);
765 return( rel_satisfied(rel) );
766 }
767 rel_set_satisfied(rel,(-res <= tolerance ));
768 return( rel_satisfied(rel) );
769 }
770 if( res > 0.0 ) {
771 if( rel_greater(rel) ) {
772 rel_set_satisfied(rel,TRUE);
773 return( rel_satisfied(rel) );
774 }
775 if( rel_less(rel) ) {
776 rel_set_satisfied(rel,FALSE);
777 return( rel_satisfied(rel) );
778 }
779 rel_set_satisfied(rel,(res <= tolerance ));
780 return( rel_satisfied(rel) );
781 }
782 rel_set_satisfied(rel,rel_equal(rel)); /* strict >0 or <0 not satisfied */
783 return( rel_satisfied(rel) );
784 }
785
786 boolean relman_calc_satisfied_scaled(struct rel_relation *rel, real64 tolerance)
787 {
788 real64 res;
789 res = rel_residual(rel);
790
791 if( res < 0.0 )
792 if( rel_less(rel) )
793 rel_set_satisfied(rel,TRUE);
794 else if( rel_greater(rel) )
795 rel_set_satisfied(rel,FALSE);
796 else
797 rel_set_satisfied(rel,(-res <= tolerance * rel_nominal(rel)));
798 else if( res > 0.0 )
799 if( rel_greater(rel) )
800 rel_set_satisfied(rel,TRUE);
801 else if( rel_less(rel) )
802 rel_set_satisfied(rel,FALSE);
803 else
804 rel_set_satisfied(rel,(res <= tolerance * rel_nominal(rel)));
805 else
806 rel_set_satisfied(rel,rel_equal(rel));
807 return( rel_satisfied(rel) );
808 }
809
810 #if REIMPLEMENT
811 real64 *relman_directly_solve_new( struct rel_relation *rel,
812 struct var_variable *solvefor, int *able, int *nsolns,
813 real64 tolerance)
814 {
815 double *value;
816 if( rel_less(rel) || rel_greater(rel) || !rel_equal(rel) ||
817 rel_extnodeinfo(rel)) {
818 *able = FALSE;
819 *nsolns = 0;
820 return(NULL);
821 }
822 else if (rel->type == e_glassbox) {
823 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
824 return value;
825 }
826 else{
827 value =
828 exprman_directly_solve(rel,rel_lhs(rel),rel_rhs(rel),
829 solvefor,able,nsolns);
830 return value;
831 }
832 }
833
834 #else /* temporary */
835 real64 *relman_directly_solve_new( struct rel_relation *rel,
836 struct var_variable *solvefor, int *able, int *nsolns, real64 tolerance)
837 {
838 double *value;
839 if( rel_less(rel) ||
840 rel_greater(rel) ||
841 !rel_equal(rel) ||
842 rel_extnodeinfo(rel)) {
843 *able = FALSE;
844 *nsolns = 0;
845 return(NULL);
846 }
847 switch (rel->type ) {
848 #if 0 /* default */
849 case e_rel_glassbox:
850 value = relman_glassbox_dsolve(rel,solvefor,able,nsolns,tolerance);
851 return value;
852 #endif
853 case e_rel_token:
854 {
855 int nvars,n;
856 const struct var_variable **vlist;
857 unsigned long vindex; /* index to the compiler */
858 nvars = rel_n_incidences(rel);
859 vlist = rel_incidence_list(rel);
860 vindex = 0;
861 for (n=0;n <nvars; n++) {
862 if (vlist[n]==solvefor) {
863 vindex = n+1; /*compiler counts from 1 */
864 break;
865 }
866 }
867 value = RelationFindRoots(IPTR(rel_instance(rel)),
868 var_lower_bound(solvefor),
869 var_upper_bound(solvefor),
870 var_nominal(solvefor),
871 tolerance,
872 &(vindex),
873 able,nsolns);
874 return value;
875 }
876 default: /* e_rel_blackbox, glassbox */
877 *able = FALSE;
878 *nsolns = 0;
879 return(NULL);
880 }
881 }
882 #endif
883
884 char *dummyrelstring(slv_system_t sys, struct rel_relation *rel, int style)
885 {
886 char *result;
887 UNUSED_PARAMETER(sys);
888 UNUSED_PARAMETER(rel);
889 UNUSED_PARAMETER(style);
890 result = ASC_NEW_ARRAY(char,80);
891 sprintf(result,"relman_make_*string_*fix not implemented yet");
892 return result;
893 }
894
895 /*
896 * converst the relations vlist into an array of int indices suitable
897 * for relation write string (compiler side indexing from 1, remember).
898 * if mORs == TRUE master indices are used, else solver.
899 */
900 static
901 void relman_cookup_indices(struct RXNameData *rd,
902 struct rel_relation *rel,int mORs)
903 {
904 int nvars,n;
905 const struct var_variable **vlist;
906
907 nvars = rel_n_incidences(rel);
908 rd->indices = rel_tmpalloc_array(1+nvars,int);
909 if (rd->indices == NULL) {
910 return;
911 }
912 vlist = rel_incidence_list(rel);
913 if (mORs) {
914 for (n = 0; n < nvars; n++) {
915 rd->indices[n+1] = var_mindex(vlist[n]);
916 }
917 } else {
918 for (n = 0; n < nvars; n++) {
919 rd->indices[n+1] = var_sindex(vlist[n]);
920 }
921 }
922 }
923
924 char *relman_make_vstring_infix(slv_system_t sys,
925 struct rel_relation *rel, int style)
926 {
927 char *sbeg;
928 struct RXNameData rd = {"x",NULL,""};
929
930 if (style) {
931 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
932 NULL,NULL,relio_ascend,NULL);
933 } else {
934 /* need to cook up output indices */
935 relman_cookup_indices(&rd,rel,1);
936 sbeg = WriteRelationString(rel_instance(rel),slv_instance(sys),
937 (WRSNameFunc)RelationVarXName,
938 &rd,relio_ascend,NULL);
939 }
940 return(sbeg);
941 }
942
943 char *relman_make_vstring_postfix(slv_system_t sys,
944 struct rel_relation *rel, int style)
945 {
946 char *sbeg;
947
948 if (style) {
949 sbeg = WriteRelationPostfixString(rel_instance(rel),slv_instance(sys));
950 } else {
951 #if REIMPLEMENT
952 left = exprman_make_xstring_postfix(rel,sys,rel_lhs(rel));
953 right = exprman_make_xstring_postfix(rel,sys,rel_rhs(rel));
954 #else
955 sbeg = ASC_NEW_ARRAY(char,60);
956 if (sbeg==NULL) return sbeg;
957 sprintf(sbeg,"relman_make_xstring_postfix not reimplemented.");
958 #endif
959 }
960 return(sbeg);
961 }

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