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

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