/[ascend]/trunk/ascend4/solver/relman.c
ViewVC logotype

Annotation of /trunk/ascend4/solver/relman.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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