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

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