/[ascend]/trunk/base/generic/compiler/relation_util.c
ViewVC logotype

Diff of /trunk/base/generic/compiler/relation_util.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 693 by johnpye, Wed Jun 21 07:00:45 2006 UTC revision 694 by johnpye, Thu Jun 22 07:54:53 2006 UTC
# Line 65  Line 65 
65  #include "rootfind.h"  #include "rootfind.h"
66  #include "func.h"  #include "func.h"
67    
68    /*------------------------------------------------------------------------------
69      DATA TYPES AND GLOBAL VARS
70    */
71    
72  int g_check_dimensions_noisy = 1;  int g_check_dimensions_noisy = 1;
73  #define GCDN g_check_dimensions_noisy  #define GCDN g_check_dimensions_noisy
74    
75  #define START 10000  /* largest power of 10 held by a short */  /** global relation pointer to avoid passing a relation recursively */
76  static struct fraction real_to_frac(double real)  static struct relation *glob_rel;
 {  
    short  num, den;  
    for( den=START; den>1 && fabs(real)*den>SHRT_MAX; den /= 10 ) ;  
    num = (short)(fabs(real)*den + 0.5);  
    if( real < 0.0 ) num = -num;  
    return( CreateFraction(num,den) );  
 }  
 #undef START  
   
77    
78  /**  /**
79      Create a static buffer to use for temporary memory storage      The following global variables are used thoughout the
80  */      functions called by RelationFindroot.
 char *tmpalloc(int nbytes)  
 /**  
  ***  Temporarily allocates a given number of bytes.  The memory need  
  ***  not be freed, but the next call to this function will reuse the  
  ***  previous allocation. Memory returned will NOT be zeroed.  
  ***  Calling with nbytes==0 will free any memory allocated.  
  **/  
 {  
   static char *ptr = NULL;  
   static int cap = 0;  
   
   if( nbytes > 0 ) {  
     if( nbytes > cap ) {  
       if( ptr ) ascfree(ptr);  
       ptr = (char *)ascmalloc(nbytes);  
       cap = nbytes;  
     }  
   }  
   else {  
     if( ptr ) ascfree(ptr);  
     ptr = NULL;  
     cap = 0;  
   }  
   return ptr;  
 }  
81    
82  /**      These should probably be located at the top of this
83      it is questionable whether this should be unified with the      file alonge with glob_rel. [OK, let it be so then. -- JP]
     ArgsForToken function in relation.c  
84  */  */
85  int ArgsForRealToken(enum Expr_enum type)  static unsigned long glob_varnum;
86  {  static int glob_done;
    switch(type) {  
    case e_zero:  
    case e_int:  
    case e_real:  
    case e_var:  
      return(0);  
   
    case e_func:  
    case e_uminus:  
       return(1);  
   
    case e_plus:  
    case e_minus:  
    case e_times:  
    case e_divide:  
    case e_power:  
    case e_ipower:  
       return(2);  
87    
88     default:  /* some data structurs...*/
       FPRINTF(ASCERR,"ArgsForRealToken: Unknown relation term type.\n");  
       return(0);  
    }  
 }  
89    
90  struct dimnode {  struct dimnode {
91     dim_type d;     dim_type d;
# Line 149  struct dimnode { Line 95  struct dimnode {
95     struct fraction power;     struct fraction power;
96  };  };
97    
98  static int IsZero(struct dimnode *node)  struct ds_soln_list {
99  {     int length,capacity;
100    if( node->type==e_zero || (node->type==e_real && node->real_const==0.0) )     double *soln;
101      return TRUE;  };
102    return FALSE;  
103  }  /*
104        Define the following if you want ASCEND to panic when it hits a
105        relation error in this file. This will help with debugging (GDB).
106    
107        Comment it out for ERROR_REPORTER to be used instead.
108    */
109    #define RELUTIL_CHECK_ABORT
110    
111    /*------------------------------------------------------------------------------
112      forward declarations
113    */
114    
115    static struct fraction real_to_frac(double real);
116    char *tmpalloc(int nbytes);
117    int ArgsForRealToken(enum Expr_enum type);
118    static int IsZero(struct dimnode *node);
119    
120    /* bunch of support functions for RelationFindRoots */
121    static double RootFind(struct relation *rel, double *lower_bound, double *upper_bound,
122        double *nominal,double *tolerance,unsigned long varnum,int *status);
123    static int CalcResidGivenValue(int *mode, int *m, unsigned long *varnum,double *val, double *u, double *f, double *g);
124    int RelationInvertTokenTop(struct ds_soln_list *soln_list);
125    int RelationInvertToken(struct relation_term **term,struct ds_soln_list *soln_list,enum safe_err *not_safe);
126    static void SetUpInvertTokenTop(struct relation_term **invert_side,double *value);
127    static int SetUpInvertToken(struct relation_term *term,struct relation_term **invert_side,double *value);
128    static int SearchEval_Branch(struct relation_term *term);
129    static void InsertBranchResult(struct relation_term *term, double value);
130    static void remove_soln( struct ds_soln_list *sl, int ndx);
131    static void append_soln( struct ds_soln_list *sl, double soln);
132    static struct relation *RelationTmpTokenCopy(CONST struct relation *src);
133    static int RelationTmpCopySide(union RelationTermUnion *old,unsigned long len,union RelationTermUnion *arr);
134    static struct relation *RelationCreateTmp(unsigned long lhslen, unsigned long rhslen, enum Expr_enum relop);
135    
136    /* the following appear only to be used locally, so I've made them static  -- JP */
137    
138    static int RelationCalcDerivative(struct Instance *i, unsigned long index, double *grad);
139    /**<
140     *  This calculates the derivative of the relation df/dx (f = lhs-rhs)
141     *  where x is the INDEX-th entry in the relation's var list.
142     *  The var list is a gl_list_t indexed from 1 to length.
143     *  Non-zero return value implies a problem.<br><br>
144     *
145     *  Notes: This function is a possible source of floating point
146     *         exceptions and should not be used during compilation.
147     */
148    
149    static enum safe_err
150    RelationCalcDerivativeSafe(struct Instance *i, unsigned long index, double *grad);
151    /**<
152     *  Calculates the derivative safely.
153     *  Non-zero return value implies a problem.
154     */
155    
156    #ifndef NDEBUG
157    static int  relutil_check_inst_and_res(struct Instance *i, double *res);
158    #endif
159    
160    #ifndef NDEBUG
161    # define CHECK_INST_RES(i,res,retval) if(!relutil_check_inst_and_res(i,res)){return retval;}
162    #else
163    # define CHECK_INST_RES(i,res,retval) ((void)0)
164    #endif
165    
166    /*------------------------------------------------------------------------------
167      SOME STUFF TO DO WITH DIMENSIONS
168    */
169    
170  /*  /*
171      @TODO what this does needs to be documented here      @TODO what this does needs to be documented here
172  */  */
173  static void apply_term_dimensions(struct relation *rel,  static void apply_term_dimensions(struct relation *rel,
174                                    struct relation_term *rt,          struct relation_term *rt,
175                                    struct dimnode *first,          struct dimnode *first,
176                                    struct dimnode *second,          struct dimnode *second,
177                                    int *con,          int *con,
178                                    int *wild)          int *wild
179  {  ){
180     enum Expr_enum type;     enum Expr_enum type;
181    
182     switch(type=RelationTermType(rt)) {     switch(type=RelationTermType(rt)) {
# Line 467  static void apply_term_dimensions(struct Line 478  static void apply_term_dimensions(struct
478     }     }
479  }  }
480    
481    /**
482        @TODO what this does needs to be documented here
483    */
484  int RelationCheckDimensions(struct relation *rel, dim_type *dimens)  int RelationCheckDimensions(struct relation *rel, dim_type *dimens)
485  {  {
486    struct dimnode *stack, *sp;    struct dimnode *stack, *sp;
# Line 567  int RelationCheckDimensions(struct relat Line 581  int RelationCheckDimensions(struct relat
581    CALCULATION FUNCTIONS    CALCULATION FUNCTIONS
582  */  */
583    
 /** global relation pointer to avoid passing a relation recursively */  
 static struct relation *glob_rel;  
   
584  /** @NOTE ANY function calling RelationBranchEvaluator should set  /** @NOTE ANY function calling RelationBranchEvaluator should set
585     glob_rel to point at the relation being evaluated.  The calling     glob_rel to point at the relation being evaluated.  The calling
586     function should also set glob_rel = NULL when it is done.     function should also set glob_rel = NULL when it is done.
# Line 2057  double CalcRelationNominal(struct Instan Line 2068  double CalcRelationNominal(struct Instan
2068        return temp;        return temp;
2069      }      }
2070    }    }
2071    if (reltype == e_blackbox)    if (reltype == e_blackbox){
2072    {      ERROR_REPORTER_HERE(ASC_PROG_ERR,"blackbox not implemented yet (assuming 1.0) (%s)",__FUNCTION__);
2073      FPRINTF(ASCERR, "error in CalcRelationNominal:\n");    }else if (reltype == e_glassbox){
2074      FPRINTF(ASCERR, "blackbox not implemented yet. Assuming 1.0\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (assuming 1.0) (%s)",__FUNCTION__);
2075    }    }else if (reltype == e_opcode){
2076    if (reltype == e_glassbox)      ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
   {  
     FPRINTF(ASCERR, "error in CalcRelationNominal:\n");  
     FPRINTF(ASCERR, "glassbox not implemented yet. Assuming 1.0\n");  
   }  
   if (reltype == e_opcode)  
   {  
     FPRINTF(ASCERR, "error in CalcRelationNominal:\n");  
     FPRINTF(ASCERR, "opcode not supported.\n");  
2077    }    }
2078    glob_rel = NULL;    glob_rel = NULL;
2079    return (double)1;    return (double)1;
# Line 2157  RelationCalcResidualPostfixSafe(struct I Line 2160  RelationCalcResidualPostfixSafe(struct I
2160    
2161      /* CONSOLE_DEBUG("..."); */      /* CONSOLE_DEBUG("..."); */
2162    
2163  #ifndef NDEBUG      CHECK_INST_RES(i,res,1);
2164      if( i == NULL ) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance");  
         return safe_problem;  
     }  
     if (res == NULL){  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relationptr");  
         return safe_problem;  
     }  
     if( InstanceKind(i) != REL_INST ) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation");  
         return safe_problem;  
     }  
 #endif  
2165      r = (struct relation *)GetInstanceRelation(i, &reltype);      r = (struct relation *)GetInstanceRelation(i, &reltype);
2166    
2167      /* CONSOLE_DEBUG("..."); */      /* CONSOLE_DEBUG("..."); */
# Line 2204  RelationCalcResidualPostfixSafe(struct I Line 2195  RelationCalcResidualPostfixSafe(struct I
2195              safe_error_to_stderr(&status);              safe_error_to_stderr(&status);
2196              break;              break;
2197          case e_blackbox:          case e_blackbox:
2198                ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Attempting to evaluate blackbox...");
2199              if ( RelationCalcResidualPostfix(i,res) != 0) {              if ( RelationCalcResidualPostfix(i,res) != 0) {
2200              status = safe_problem;              status = safe_problem;
2201              safe_error_to_stderr(&status);              safe_error_to_stderr(&status);
# Line 2235  RelationCalcResidualPostfix(struct Insta Line 2227  RelationCalcResidualPostfix(struct Insta
2227    struct relation *r;    struct relation *r;
2228    enum Expr_enum reltype;    enum Expr_enum reltype;
2229    
2230  #ifndef NDEBUG    CHECK_INST_RES(i,res,1);
2231    if( i == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualPostfix: null instance\n");  
     return 1;  
   }  
   if (res == NULL){  
     FPRINTF(ASCERR,"error in RelationCalcResidualPostfix: null relationptr\n");  
     return 1;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualPostfix: not relation\n");  
     return 1;  
   }  
 #endif  
2232    r = (struct relation *)GetInstanceRelation(i, &reltype);    r = (struct relation *)GetInstanceRelation(i, &reltype);
2233    if( r == NULL ) {    if( r == NULL ) {
2234      FPRINTF(ASCERR, "error in RelationCalcResidualPostfix: null relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2235      return 1;      return 1;
2236    }    }
2237    if( reltype == e_token ) {    if( reltype == e_token ) {
# Line 2272  RelationCalcResidualPostfix(struct Insta Line 2252  RelationCalcResidualPostfix(struct Insta
2252      }      }
2253      return 0;      return 0;
2254    } else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    } else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2255      if (reltype == e_blackbox)  
2256      {      if (reltype == e_blackbox){
2257        /* FIXME */        /* FIXME */
2258      /* note: blackbox equations support the form        /* note: blackbox equations support the form
2259          output[i] = f(input[j] for all j) foreach i           output[i] = f(input[j] for all j) foreach i
2260      thus the residual is           thus the residual is ... (?)
2261      */        */
2262          ERROR_REPORTER_HERE(ASC_PROG_ERR,"blackbox not implemented yet (%s)",__FUNCTION__);
2263        FPRINTF(ASCERR, "error in RelationCalcResidualPostfix:\n");      }else if (reltype == e_glassbox){
2264        FPRINTF(ASCERR, "blackbox not implemented yet.\n");        ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2265      }      }else if (reltype == e_opcode)    {
2266      if (reltype == e_glassbox)        ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
     {  
       FPRINTF(ASCERR, "error in RelationCalcResidualPostfix:\n");  
       FPRINTF(ASCERR, "glassbox not implemented yet.\n");  
     }  
     if (reltype == e_opcode)  
     {  
       FPRINTF(ASCERR, "error in RelationCalcResidualPostfix:\n");  
       FPRINTF(ASCERR, "opcode not supported.\n");  
2267      }      }
2268    
2269      return 1;      return 1;
2270    } else {    } else {
2271      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__,"reached end of routine");
               "error in RelationCalcResidualPostfix:\n"  
               "reached end of routine\n");  
2272      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2273    }    }
2274  }  }
# Line 2310  int RelationCalcExceptionsInfix(struct I Line 2281  int RelationCalcExceptionsInfix(struct I
2281    int old_errno;    int old_errno;
2282    
2283    glob_rel = NULL;    glob_rel = NULL;
2284  #ifndef NDEBUG  
2285    if( i == NULL ) {    CHECK_INST_RES(i,&res,-1);
2286      FPRINTF(ASCERR, "error in RelationCalcExceptionsInfix: NULL instance\n");  
     return -1;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcExceptionsInfix: not relation\n");  
     return -1;  
   }  
 #endif  
2287    glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);    glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);
2288    if( glob_rel == NULL ) {    if( glob_rel == NULL ) {
2289      FPRINTF(ASCERR, "error in RelationCalcExceptionsInfix: NULL relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation");
2290      return -1;      return -1;
2291    }    }
2292    if( reltype == e_token ) {    if( reltype == e_token ) {
# Line 2359  int RelationCalcExceptionsInfix(struct I Line 2323  int RelationCalcExceptionsInfix(struct I
2323      }      }
2324      glob_rel = NULL;      glob_rel = NULL;
2325      return result;      return result;
2326    } else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2327      FPRINTF(ASCERR, "error in RelationCalcExceptionsInfix:\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"relation type not implemented (%s)",__FUNCTION__);
     FPRINTF(ASCERR, "reltype not implemented yet\n");  
2328      glob_rel = NULL;      glob_rel = NULL;
2329      return -1;      return -1;
2330    } else {    }else{
2331      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__,"reached end of routine\n");
               "error in RelationCalcExceptionsInfix:\n"  
               "reached end of routine\n");  
2332      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2333    }    }
2334  }  }
2335    
2336    
2337  int RelationCalcResidualInfix(struct Instance *i, double *res)  int RelationCalcResidualInfix(struct Instance *i, double *res)
2338  {  {
2339    enum Expr_enum reltype;    enum Expr_enum reltype;
2340    glob_rel = NULL;    glob_rel = NULL;
2341    
2342  #ifndef NDEBUG    CHECK_INST_RES(i,res,1);
2343    if( i == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualInfix: NULL instance\n");  
     return 1;  
   }  
   if (res == NULL){  
     FPRINTF(ASCERR,"error in RelationCalcResidualInfix: NULL residual ptr\n");  
     return 1;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualInfix: not relation\n");  
     return 1;  
   }  
 #endif  
2344    glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);    glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);
2345    if( glob_rel == NULL ) {    if( glob_rel == NULL ) {
2346      FPRINTF(ASCERR, "error in RelationCalcResidualInfix: NULL relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation\n");
2347      return 1;      return 1;
2348    }    }
2349    if( reltype == e_token ) {    if( reltype == e_token ) {
# Line 2407  int RelationCalcResidualInfix(struct Ins Line 2357  int RelationCalcResidualInfix(struct Ins
2357      }      }
2358      glob_rel = NULL;      glob_rel = NULL;
2359      return 0;      return 0;
2360    } else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2361      FPRINTF(ASCERR, "error in RelationCalcResidualInfix:\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
     FPRINTF(ASCERR, "reltype not implemented yet\n");  
2362      glob_rel = NULL;      glob_rel = NULL;
2363      return 1;      return 1;
2364    } else {    }else{
2365      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__,"reached end of routine\n");
               "error in RelationCalcResidualInfix:\n"  
               "reached end of routine\n");  
2366      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2367    }    }
2368  }  }
# Line 2431  RelationCalcResidualPostfix2(struct Inst Line 2378  RelationCalcResidualPostfix2(struct Inst
2378    struct relation *r;    struct relation *r;
2379    enum Expr_enum reltype;    enum Expr_enum reltype;
2380    
2381  #ifndef NDEBUG    CHECK_INST_RES(i,res,1);
2382    if( i == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualPostfix2: null instance\n");  
     return 1;  
   }  
   if( res == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualPostfix2: %s\n",  
             "null relation ptr");  
     return 1;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidualPostfix2: not relation\n");  
     return 1;  
   }  
 #endif  
2383    r = (struct relation *)GetInstanceRelation(i, &reltype);    r = (struct relation *)GetInstanceRelation(i, &reltype);
2384    if( r == NULL ) {    if( r == NULL ) {
2385      FPRINTF(ASCERR, "error in RelationCalcResidualPostfix2: null relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2386      return 1;      return 1;
2387    }    }
2388    
2389    if( reltype == e_token ) {    if( reltype == e_token ){
2390      *res = RelationEvaluateResidualPostfix(r);      *res = RelationEvaluateResidualPostfix(r);
2391      return 0;      return 0;
2392    } else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2393      FPRINTF(ASCERR, "error in RelationCalcResidualPostfix2:\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
     FPRINTF(ASCERR, "reltype not implemented yet\n");  
2394      return 1;      return 1;
2395    } else {    }else{
2396      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__,"reached end of routine\n");
               "error in RelationCalcResidualPostfix2:\n"  
               "reached end of routine\n");  
2397      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2398    }    }
2399  }  }
# Line 2502  RelationCalcResidGrad(struct Instance *i Line 2433  RelationCalcResidGrad(struct Instance *i
2433    struct relation *r;    struct relation *r;
2434    enum Expr_enum reltype;    enum Expr_enum reltype;
2435    
2436  #ifndef NDEBUG    CHECK_INST_RES(i,residual,1);
2437    if( i == NULL ) {    CHECK_INST_RES(i,residual,1);
2438      FPRINTF(ASCERR, "error in RelationCalcResidGrad: null instance\n");  
     return 1;  
   }  
   if( residual == NULL || gradient == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidGrad: passed a null pointer\n");  
     return 1;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcResidGrad: not relation\n");  
     return 1;  
   }  
 #endif  
2439    r = (struct relation *)GetInstanceRelation(i, &reltype);    r = (struct relation *)GetInstanceRelation(i, &reltype);
2440    if( r == NULL ) {    if( r == NULL ) {
2441      FPRINTF(ASCERR, "error in RelationCalcResidGrad: null relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2442      return 1;      return 1;
2443    }    }
2444    
2445    if( reltype == e_token ) {    if(reltype == e_token ){
2446      return RelationEvaluateResidualGradient(r, residual, gradient);      return RelationEvaluateResidualGradient(r, residual, gradient);
2447    }  
2448    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    }else if(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2449      FPRINTF(ASCERR, "error in RelationCalcResidGrad %s\n",      ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
             "reltype not implemented");  
2450      return 1;      return 1;
2451    }  
2452    else {    }else{
2453      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__, "reached end of routine");
               "error in RelationCalcResidGrad:\n"  
               "reached end of routine");  
2454      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2455    }    }
2456  }  }
# Line 2550  RelationCalcResidGradSafe(struct Instanc Line 2467  RelationCalcResidGradSafe(struct Instanc
2467    
2468  #ifndef NDEBUG  #ifndef NDEBUG
2469    if( i == NULL ) {    if( i == NULL ) {
2470      FPRINTF(ASCERR, "error in RelationCalcResidGradSafe: null instance\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance\n");
2471      not_safe = safe_problem;      not_safe = safe_problem;
2472      return not_safe;      return not_safe;
2473    }    }
2474    if( residual == NULL || gradient == NULL ) {    if( residual == NULL || gradient == NULL ) {
2475      FPRINTF(ASCERR,      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null pointer\n");
             "error in RelationCalcResidGradSafe: passed a null pointer\n");  
2476      not_safe = safe_problem;      not_safe = safe_problem;
2477      return not_safe;      return not_safe;
2478    }    }
2479    if( InstanceKind(i) != REL_INST ) {    if( InstanceKind(i) != REL_INST ) {
2480      FPRINTF(ASCERR, "error in RelationCalcResidGradSafe: not relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation\n");
2481      not_safe = safe_problem;      not_safe = safe_problem;
2482      return not_safe;      return not_safe;
2483    }    }
2484  #endif  #endif
2485    r = (struct relation *)GetInstanceRelation(i, &reltype);    r = (struct relation *)GetInstanceRelation(i, &reltype);
2486    if( r == NULL ) {    if( r == NULL ) {
2487      FPRINTF(ASCERR, "error in RelationCalcResidGradSafe: null relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2488      not_safe = safe_problem;      not_safe = safe_problem;
2489      return not_safe;      return not_safe;
2490    }    }
# Line 2579  RelationCalcResidGradSafe(struct Instanc Line 2495  RelationCalcResidGradSafe(struct Instanc
2495      return not_safe;      return not_safe;
2496    }    }
2497    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2498      if (reltype == e_blackbox)      if (reltype == e_blackbox){
2499      {        ERROR_REPORTER_HERE(ASC_PROG_ERR,"blackbox not implemented yet (%s)",__FUNCTION__);
2500        FPRINTF(ASCERR, "error in RelationCalcResidGradSafe:\n");      }else if (reltype == e_glassbox){
2501        FPRINTF(ASCERR, "blackbox not implemented yet.\n");        ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2502      }      }else if (reltype == e_opcode){
2503      if (reltype == e_glassbox)        ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
     {  
       FPRINTF(ASCERR, "error in RelationCalcResidGradSafe:\n");  
       FPRINTF(ASCERR, "glassbox not implemented yet.\n");  
     }  
     if (reltype == e_opcode)  
     {  
       FPRINTF(ASCERR, "error in RelationCalcResidGradSafe:\n");  
       FPRINTF(ASCERR, "opcode not supported.\n");  
2504      }      }
2505      not_safe = safe_problem;      not_safe = safe_problem;
2506      return not_safe;      return not_safe;
# Line 2609  RelationCalcResidGradSafe(struct Instanc Line 2517  RelationCalcResidGradSafe(struct Instanc
2517  /*  /*
2518      calculate the derivative with respect to a single variable      calculate the derivative with respect to a single variable
2519      whose index is index, where 1<=index<=NumberVariables(r)      whose index is index, where 1<=index<=NumberVariables(r)
2520    
2521        @TODO this appears only to be used in PrintGradients
2522  */  */
2523  int  int
2524  RelationCalcDerivative(struct Instance *i,  RelationCalcDerivative(struct Instance *i,
# Line 2618  RelationCalcDerivative(struct Instance * Line 2528  RelationCalcDerivative(struct Instance *
2528    struct relation *r;    struct relation *r;
2529    enum Expr_enum reltype;    enum Expr_enum reltype;
2530    
2531  #ifndef NDEBUG    CHECK_INST_RES(i,gradient,1);
2532    if( i == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcDerivative: null instance\n");  
     return 1;  
   }  
   if( gradient == NULL ) {  
     FPRINTF(ASCERR, "error in RelationCalcDerivative: passed null pointer\n");  
     return 1;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcDerivative: not relation\n");  
     return 1;  
   }  
 #endif  
2533    r = (struct relation *)GetInstanceRelation(i, &reltype);    r = (struct relation *)GetInstanceRelation(i, &reltype);
2534    if( r == NULL ) {    if( r == NULL ) {
2535      FPRINTF(ASCERR, "error in RelationCalcDerivative: null relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2536      return 1;      return 1;
2537    }    }
2538    if( (index < 1) || (index > NumberVariables(r)) ) {    if( (index < 1) || (index > NumberVariables(r)) ) {
2539      FPRINTF(ASCERR, "error in RelationCalcDerivative: index out of bounds\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"index out of bounds\n");
2540      return 1;      return 1;
2541    }    }
2542    
# Line 2647  RelationCalcDerivative(struct Instance * Line 2545  RelationCalcDerivative(struct Instance *
2545      return 0;      return 0;
2546    }    }
2547    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2548      FPRINTF(ASCERR, "error in RelationCalcDerivative %s\n",      ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not supported (%s)",__FUNCTION__);
             "reltype not implemented");  
2549      return 1;      return 1;
2550    }    }
2551    else {    else {
2552      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__,"reached end of routine");
               "error in RelationCalcDerivative: \n"  
               "reached end of routine");  
2553      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2554    }    }
2555  }  }
# Line 2669  RelationCalcDerivativeSafe(struct Instan Line 2564  RelationCalcDerivativeSafe(struct Instan
2564    enum safe_err not_safe = safe_ok;    enum safe_err not_safe = safe_ok;
2565    
2566  #ifndef NDEBUG  #ifndef NDEBUG
2567    if( i == NULL ) {    if(!relutil_check_inst_and_res(i,gradient)){
     FPRINTF(ASCERR, "error in RelationCalcDerivativeSafe: null instance\n");  
2568      not_safe = safe_problem;      not_safe = safe_problem;
2569      return not_safe;      return not_safe;
   }  
   if( gradient == NULL ) {  
     FPRINTF(ASCERR,  
             "error in RelationCalcDerivativeSafe: passed null pointer\n");  
     not_safe = safe_problem;  
     return not_safe;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationCalcDerivativeSafe: not relation\n");  
     not_safe = safe_problem;  
     return not_safe;  
2570    }    }
2571  #endif  #endif
2572    r = (struct relation *)GetInstanceRelation(i, &reltype);    r = (struct relation *)GetInstanceRelation(i, &reltype);
2573    if( r == NULL ) {    if( r == NULL ) {
2574      FPRINTF(ASCERR, "error in RelationCalcDerivativeSafe: null relation\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2575      not_safe = safe_problem;      not_safe = safe_problem;
2576      return not_safe;      return not_safe;
2577    }    }
2578    if( (index < 1) || (index > NumberVariables(r)) ) {    if( (index < 1) || (index > NumberVariables(r)) ) {
2579      FPRINTF(ASCERR,      ERROR_REPORTER_HERE(ASC_PROG_ERR,"index out of bounds\n");
             "error in RelationCalcDerivativeSafe: index out of bounds\n");  
2580      not_safe = safe_problem;      not_safe = safe_problem;
2581      return not_safe;      return not_safe;
2582    }    }
# Line 2704  RelationCalcDerivativeSafe(struct Instan Line 2586  RelationCalcDerivativeSafe(struct Instan
2586      return not_safe;      return not_safe;
2587    }    }
2588    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {    else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2589      FPRINTF(ASCERR, "error in RelationCalcDerivativeSafe %s\n",      ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not supported (%s)",__FUNCTION__);
             "reltype not implemented");  
2590      not_safe = safe_problem;      not_safe = safe_problem;
2591      return not_safe;      return not_safe;
2592    }    }
2593    else {    else {
2594      Asc_Panic(2, NULL,      Asc_Panic(2, __FUNCTION__, "reached end of routine");
               "error in RelationCalcDerivativeSafe:\n"  
               "reached end of routine");  
2595      exit(2);/* Needed to keep gcc from whining */      exit(2);/* Needed to keep gcc from whining */
2596    }    }
2597  }  }
# Line 2913  void PrintRelationResiduals(struct Insta Line 2792  void PrintRelationResiduals(struct Insta
2792    VisitInstanceTree(i,PrintResidual, 0, 0);    VisitInstanceTree(i,PrintResidual, 0, 0);
2793  }  }
2794    
2795    /*==============================================================================
2796      'RELATIONFINDROOTS' AND SUPPORT FUNCTIONS
2797    
2798  /**      The following functions support RelationFindRoots which
2799   *** The following functions support RelationFindRoots which      is the compiler implementation of our old DirectSolve
2800   *** is the compiler implementation of our old DirectSolve      function.
  *** function.  These functions can be catagorized as follows:  
  *** Memory Management and Copying functions:  
  ***      RelationCreateTmp, RelationTmpCopySide,  
  ***      RelationTmpTokenCopy, append_soln, remove_soln  
  *** Direct Solve Functions:  
  ***      InsertBranchResult, SearchEval_Branch, SetUpInvertToken,  
  ***      SetUpInvertTokenTop, RelationInvertToken, RelationInvertTokenTop  
  *** Rootfinding Functions:  
  ***      CalcResidGivenValue, RootFind  
  *** Eternal Function:  
  ***      RelationFindRoots  
  **/  
2801    
2802  /*------------------------------------------------------------------------------      I suspect that a lot of this stuff is deprecated -- JP
   MEMORY MANAGEMENT AND COPYING FUNCTIONS  
2803  */  */
2804    
2805  /*  double *RelationFindRoots(struct Instance *i,
2806   * RelationCreateTmp creates a struct relation of type e_token          double lower_bound, double upper_bound,
2807   * and passes back a pointer to the relation.  The lengths of          double nominal,
2808   * the right and left sides (lhslen and rhslen) of the relation          double tolerance,
2809   * are supplied by the calling function.          unsigned long *varnum,
2810   * User is responsible for setting RTOKEN(return).*_len.          int *able,
2811   * Basically, all this does is manage memory nicely.          int *nsolns
2812   *  ){
2813   * IF called with all 0/NULL, frees internal recycles.    struct relation *rel;
2814   */    double sideval;
2815  static    enum Expr_enum reltype;
2816  struct relation *RelationCreateTmp(unsigned long lhslen, unsigned long rhslen,    static struct ds_soln_list soln_list = {0,0,NULL};
2817                                     enum Expr_enum relop)    CONST struct gl_list_t *list;
 {  
   static struct relation *rel=NULL;  
   static unsigned long lhscap=0, rhscap=0;  
2818    
2819    /* check for recycle clear and free things if needed. */    /* check for recycle shutdown */
2820    if (lhslen==0 && rhslen == 0 && relop == e_nop) {    if (i==NULL && varnum == NULL && able == NULL && nsolns == NULL) {
2821      if (rel != NULL) {      if (soln_list.soln != NULL) {
2822        if (rel->share != NULL) {        ascfree(soln_list.soln);
2823          if (RTOKEN(rel).lhs!=NULL) {        soln_list.soln = NULL;
2824            ascfree(RTOKEN(rel).lhs);        soln_list.length = soln_list.capacity = 0;
         }  
         if (RTOKEN(rel).rhs!=NULL)  {  
           ascfree(RTOKEN(rel).rhs);  
         }  
         ascfree(rel->share);  
       }  
       ascfree(rel);  
       rel = NULL;  
2825      }      }
2826      lhscap = rhscap = 0;      RootFind(NULL,NULL,NULL,NULL,NULL,0L,NULL); /*clear brent recycle */
2827        RelationCreateTmp(0,0,e_nop); /* clear tmprelation recycle */
2828      return NULL;      return NULL;
2829    }    }
2830    if (rel == NULL) {    /* check assertions */
2831      rel = CreateRelationStructure(relop,crs_NEWUNION);  #ifndef NDEBUG
2832    }    if( i == NULL ) {
2833    if (lhscap < lhslen) {      FPRINTF(ASCERR, "error in RelationFindRoot: NULL instance\n");
2834      lhscap = lhslen;      glob_rel = NULL;
2835      if ( RTOKEN(rel).lhs != NULL) {      return NULL;
2836        ascfree(RTOKEN(rel).lhs);    }
2837      if (able == NULL){
2838        FPRINTF(ASCERR,"error in RelationFindRoot: NULL able ptr\n");
2839        glob_rel = NULL;
2840        return NULL;
2841      }
2842      if (varnum == NULL){
2843        FPRINTF(ASCERR,"error in RelationFindRoot: NULL varnum\n");
2844        glob_rel = NULL;
2845        return NULL;
2846      }
2847      if( InstanceKind(i) != REL_INST ) {
2848        FPRINTF(ASCERR, "error in RelationFindRoot: not relation\n");
2849        glob_rel = NULL;
2850        return NULL;
2851      }
2852    #endif
2853    
2854      *able = FALSE;
2855      *nsolns = -1;     /* nsolns will be -1 for a very unhappy root-finder */
2856      glob_rel = NULL;
2857      glob_done = 0;
2858      soln_list.length = 0; /* reset len to 0. if NULL to start, append mallocs */
2859      append_soln(&soln_list,0.0);
2860      rel = (struct relation *)GetInstanceRelation(i, &reltype);
2861      if( rel == NULL ) {
2862        FPRINTF(ASCERR, "error in RelationFindRoot: NULL relation\n");
2863        glob_rel = NULL; return NULL;
2864      }
2865      /* here we should switch and handle all types. at present we don't
2866       * handle anything except e_token
2867       */
2868      if( reltype != e_token ) {
2869        FPRINTF(ASCERR, "error in RelationFindRoot: non-token relation\n");
2870        glob_rel = NULL;
2871        return NULL;
2872      }
2873    
2874      if (RelationRelop(rel) == e_equal){
2875        glob_rel = RelationTmpTokenCopy(rel);
2876        assert(glob_rel!=NULL);
2877        glob_done = 0;
2878        list = RelationVarList(glob_rel);
2879        if( *varnum >= 1 && *varnum <= gl_length(list)){
2880          glob_done = 1;
2881        }
2882        if (!glob_done) {
2883          FPRINTF(ASCERR, "error in FindRoot: var not found\n");
2884          glob_rel = NULL;
2885          return NULL;
2886        }
2887    
2888        glob_varnum = *varnum;
2889        glob_done = 0;
2890        assert(Infix_LhsSide(glob_rel) != NULL);
2891        /* In the following if statements we look for the target variable
2892         * to the left and right, evaluating all branches without the
2893         * target.
2894         */
2895        if (SearchEval_Branch(Infix_LhsSide(glob_rel)) < 1) {
2896          sideval = RelationBranchEvaluator(Infix_LhsSide(glob_rel));
2897          if (finite(sideval)) {
2898            InsertBranchResult(Infix_LhsSide(glob_rel),sideval);
2899          } else {
2900            FPRINTF(ASCERR,"Inequality in RelationFindRoots. Infinite RHS.\n");
2901            glob_rel = NULL;
2902            return NULL;
2903          }
2904        }
2905        assert(Infix_RhsSide(glob_rel) != NULL);
2906        if (SearchEval_Branch(Infix_RhsSide(glob_rel)) < 1) {
2907            sideval = RelationBranchEvaluator(Infix_RhsSide(glob_rel));
2908            if (finite(sideval)) {
2909              InsertBranchResult(Infix_RhsSide(glob_rel),sideval);
2910            } else {
2911              FPRINTF(ASCERR,"Inequality in RelationFindRoots. Infinite LHS.\n");
2912              glob_rel = NULL;
2913              return NULL;
2914            }
2915        }
2916        if (glob_done < 1) {
2917          /* RelationInvertToken never found variable */
2918          glob_done = 0;
2919          *able = FALSE;
2920          return soln_list.soln;
2921        }
2922        if (glob_done == 1) {
2923          /* set to 0 so while loop in RelationInvertToken will work */
2924          glob_done = 0;
2925          glob_done = RelationInvertTokenTop(&(soln_list));
2926        }
2927        if (glob_done == 1) { /* if still one, token inversions successful */
2928          glob_done = 0;
2929          *nsolns= soln_list.length;
2930          *able = TRUE;
2931          return soln_list.soln;
2932        }
2933        /* CALL ITERATIVE SOLVER */
2934        *soln_list.soln = RootFind(glob_rel,&(lower_bound),
2935                           &(upper_bound),&(nominal),
2936                           &(tolerance),
2937                           glob_varnum,able);
2938    
2939        glob_done = 0;
2940        if(*able == 0) { /* Root-Find returns 0 for success*/
2941          *nsolns = 1;
2942          *able = TRUE;
2943        } else {
2944          *able = FALSE;
2945        }
2946        return soln_list.soln;
2947    
2948      }
2949      FPRINTF(ASCERR,"Inequality in RelationFindRoots. can't find roots.\n");
2950      *able = FALSE;
2951      return soln_list.soln;
2952    }
2953    
2954    /*------------------------------------------------------------------------------
2955      MEMORY MANAGEMENT AND COPYING FUNCTIONS
2956      to support the RelationFindRoots function.
2957    */
2958    
2959    /**
2960        @see RelationFindRoots
2961    
2962        Create a struct relation of type e_token
2963        and passes back a pointer to the relation.
2964        
2965        The lengths of
2966        the right and left sides (lhslen and rhslen) of the relation
2967        are supplied by the calling function.
2968        
2969        User is responsible for setting RTOKEN(return).*_len.
2970        
2971        Basically, all this does is manage memory nicely.
2972    
2973        IF called with all 0/NULL, frees internal recycles.
2974    */
2975    static struct relation *RelationCreateTmp(
2976            unsigned long lhslen, unsigned long rhslen,
2977            enum Expr_enum relop
2978    ){
2979      static struct relation *rel=NULL;
2980      static unsigned long lhscap=0, rhscap=0;
2981    
2982      /* check for recycle clear and free things if needed. */
2983      if (lhslen==0 && rhslen == 0 && relop == e_nop) {
2984        if (rel != NULL) {
2985          if (rel->share != NULL) {
2986            if (RTOKEN(rel).lhs!=NULL) {
2987              ascfree(RTOKEN(rel).lhs);
2988            }
2989            if (RTOKEN(rel).rhs!=NULL)  {
2990              ascfree(RTOKEN(rel).rhs);
2991            }
2992            ascfree(rel->share);
2993          }
2994          ascfree(rel);
2995          rel = NULL;
2996        }
2997        lhscap = rhscap = 0;
2998        return NULL;
2999      }
3000      if (rel == NULL) {
3001        rel = CreateRelationStructure(relop,crs_NEWUNION);
3002      }
3003      if (lhscap < lhslen) {
3004        lhscap = lhslen;
3005        if ( RTOKEN(rel).lhs != NULL) {
3006          ascfree(RTOKEN(rel).lhs);
3007      }      }
3008      RTOKEN(rel).lhs = (union RelationTermUnion *)      RTOKEN(rel).lhs = (union RelationTermUnion *)
3009          ascmalloc(lhscap*sizeof(union RelationTermUnion));          ascmalloc(lhscap*sizeof(union RelationTermUnion));
# Line 2992  struct relation *RelationCreateTmp(unsig Line 3020  struct relation *RelationCreateTmp(unsig
3020  }  }
3021    
3022  /**  /**
3023   *** The following global variables are used thoughout the      @see RelationFindRoots
  *** functions called by RelationFindroot.  
  *** These should probably be located at the top of this  
  *** file alonge with glob_rel.  
  **/  
 static unsigned long glob_varnum;  
 static int glob_done;  
3024    
3025  /**      Full-blown relation copy (not copy by reference)
3026   *** The following is documentation from the old exprman  
3027   *** file.  RelationTmpCopySide and RelationTmpCopyToken      We can now just do a memcopy and the infix pointers
3028   *** are reimplimentations of exprman functions.      all adjust by the difference between the token
3029   **/      arrays that the gl_lists are hiding. Cool, eh?
3030  /*      
3031   * We can now just do a memcopy and the infix pointers      @NOTE if any turkey ever tries to delete an individual
3032   * all adjust by the difference between the token      token from these gl_lists AND deallocate it,
3033   * arrays that the gl_lists are hiding. Cool, eh?      they will get a severe headache. Ooo scary.
3034   * Note, if any turkey ever tries to delete an individual  
3035   * token from these gl_lists AND deallocate it,      This is a full blown copy and not copy by reference.
3036   * they will get a severe headache.      You do not need to remake the infix pointers after
3037   *      calling this function. return 0 if ok, 1 if error.
3038   * This is a full blown copy and not copy by reference.  
3039   * You do not need to remake the infix pointers after      @NOTE RelationTmpCopySide and RelationTmpCopyToken are reimplimentations
3040   * calling this function. return 0 if ok, 1 if error.      of functions from the v. old 'exprman' file.
3041   */  */
3042  static int RelationTmpCopySide(union RelationTermUnion *old,  static int RelationTmpCopySide(union RelationTermUnion *old,
3043                      unsigned long len,          unsigned long len,
3044                      union RelationTermUnion *arr)          union RelationTermUnion *arr
3045  {  ){
3046    struct relation_term *term;    struct relation_term *term;
3047    unsigned long c;    unsigned long c;
3048    long int delta;    long int delta;
# Line 3080  static int RelationTmpCopySide(union Rel Line 3102  static int RelationTmpCopySide(union Rel
3102    return 0;    return 0;
3103  }  }
3104    
3105  /*  /**
3106   * The relation returned by this function should have      @see RelationFindRoots
3107   * NO persistent pointers made to it, as it is still  
3108   * our property. The vars in the relation do not      Copy tmp token for a relation (guess -- JP)
3109   * know about these references to them, as this is  
3110   * a tmp rel.      The relation returned by this function should have
3111   */      NO persistent pointers made to it, as it is still
3112  static      our property. The vars in the relation do not
3113  struct relation *RelationTmpTokenCopy(CONST struct relation *src)      know about these references to them, as this is
3114  {      a tmp rel.
3115    
3116        @NOTE RelationTmpCopySide and RelationTmpCopyToken are reimplimentations
3117        of functions from the v. old 'exprman' file.
3118    */
3119    static struct relation *RelationTmpTokenCopy(CONST struct relation *src){
3120    struct relation *result;    struct relation *result;
3121    long int delta;    long int delta;
3122    assert(src!=NULL);    assert(src!=NULL);
# Line 3125  struct relation *RelationTmpTokenCopy(CO Line 3152  struct relation *RelationTmpTokenCopy(CO
3152    return result;    return result;
3153  }  }
3154    
 struct ds_soln_list {  
    int length,capacity;  
    double *soln;  
 };  
   
   
3155  #define alloc_array(nelts,type)   \  #define alloc_array(nelts,type)   \
3156     ((nelts) > 0 ? (type *)ascmalloc((nelts)*sizeof(type)) : NULL)     ((nelts) > 0 ? (type *)ascmalloc((nelts)*sizeof(type)) : NULL)
3157  #define copy_nums(from,too,nnums)  \  #define copy_nums(from,too,nnums)  \
3158     asc_memcpy((from),(too),(nnums)*sizeof(double))     asc_memcpy((from),(too),(nnums)*sizeof(double))
3159    
 static  
 void append_soln( struct ds_soln_list *sl, double soln)  
3160  /**  /**
3161   ***  Appends the solution onto the solution list      @see RelationFindRoots
3162   **/  
3163  {      Appends the solution onto the solution list
3164    */
3165    static void append_soln( struct ds_soln_list *sl, double soln){
3166    if( sl->length == sl->capacity ) {    if( sl->length == sl->capacity ) {
3167      int newcap;      int newcap;
3168      double *newlist;      double *newlist;
# Line 3159  void append_soln( struct ds_soln_list *s Line 3180  void append_soln( struct ds_soln_list *s
3180    sl->soln[sl->length++] = soln;    sl->soln[sl->length++] = soln;
3181  }  }
3182    
3183  static  /**
3184  void remove_soln( struct ds_soln_list *sl, int ndx)      @see RelationFindRoots
3185  /*  
3186   *  Removes solution at given index from solution list.      Removes solution at given index from solution list.
3187   */  */
3188  {  static void remove_soln( struct ds_soln_list *sl, int ndx){
3189    copy_nums((char *)(sl->soln+ndx+1),    copy_nums((char *)(sl->soln+ndx+1),
3190              (char *)(sl->soln+ndx), --(sl->length) - ndx);              (char *)(sl->soln+ndx), --(sl->length) - ndx);
3191  }  }
3192    
3193  /*------------------------------------------------------------------------------  /*------------------------------------------------------------------------------
3194    DIRECT SOLVE FUNCTIONS    DIRECT SOLVE FUNCTIONS
3195      to support the RelationFindRoots function.
3196  */  */
3197    
3198  /**  /**
3199   *** InsertBranchResult changes a relation term type to e_real and      @see RelationFindRoots
3200   *** fills the value field of this term.  In subsequent passes of  
3201   *** the RelationBranchEvaluator the term will be considered to      Change a relation term type to e_real and fill the value field of this term.
3202   *** be a leaf.  
3203        In subsequent passes of the RelationBranchEvaluator the term will be
3204        considered to be a leaf.
3205   */   */
3206  static void InsertBranchResult(struct relation_term *term, double value)  static void InsertBranchResult(struct relation_term *term, double value){
 {  
3207    assert(term!=NULL);    assert(term!=NULL);
3208    term->t = e_real;    term->t = e_real;
3209    R_TERM(term)->value = value;    R_TERM(term)->value = value;
3210  }  }
3211    
3212  /**  /**
3213   *** SearchEval_Branch simplifies branches of a relation      @see RelationFindRoots
3214   *** (the relation pointed to by glob_rel).  Only terms  
3215   *** of type e_real, e_int, e_zero, and e_var are left      Simplify branches of a relation (the relation pointed to by glob_rel).
3216   *** hanging off the operators on the path to the  
3217   *** variable (with varnum = glob_varnum) being direct      Only terms of type e_real, e_int, e_zero, and e_var are left
3218   *** solved for.      hanging off the operators on the path to the
3219   *** This may need to be changed to only leave e_reals      variable (with varnum = glob_varnum) being direct
3220   *** so that the inversion routine can make faster decisions???      solved for.
3221   *** Probably not.  
3222   ***      @TODO This may need to be changed to only leave e_reals
3223   *** Returns >= 1 if glob_varnum spotted, else 0 (or at least <1).      so that the inversion routine can make faster decisions???
3224   **/      Probably not.
3225  static int SearchEval_Branch(struct relation_term *term)      
3226  {      @return >= 1 if glob_varnum spotted, else 0 (or at least <1).
3227    */
3228    static int SearchEval_Branch(struct relation_term *term){
3229    int result = 0;    int result = 0;
3230    assert(term != NULL);    assert(term != NULL);
3231    switch(RelationTermType(term)) {    switch(RelationTermType(term)) {
# Line 3278  static int SearchEval_Branch(struct rela Line 3303  static int SearchEval_Branch(struct rela
3303  }  }
3304    
3305  /**  /**
3306   *** SetUpInvertToken selects the side of the relation which      @see RelationFindRoots
3307   *** will be inverted next.  It also fills the value which is  
3308   *** currently being inverted on.      Select the side of the relation which will be inverted next.
3309   *** This function assumes SearchEval_Branch has been called  
3310   *** previously.      Also fill the value which is currently being inverted on.
3311   **/      
3312        @NOTE This function assumes SearchEval_Branch has been called
3313        previously.
3314    */
3315  static int SetUpInvertToken(struct relation_term *term,  static int SetUpInvertToken(struct relation_term *term,
3316                   struct relation_term **invert_side,          struct relation_term **invert_side,
3317                   double *value)          double *value
3318  {  ){
3319    switch(RelationTermType(term)) {    switch(RelationTermType(term)) {
3320    case e_uminus:    case e_uminus:
3321        *invert_side = TermBinLeft(term);        *invert_side = TermBinLeft(term);
# Line 3323  static int SetUpInvertToken(struct relat Line 3351  static int SetUpInvertToken(struct relat
3351    }    }
3352  }  }
3353    
3354    /**
3355  static void SetUpInvertTokenTop(struct relation_term **invert_side,      @see RelationFindRoots
3356                   double *value)  */
3357  {  static void SetUpInvertTokenTop(
3358            struct relation_term **invert_side,
3359            double *value
3360    ){
3361    switch(RelationTermType(Infix_RhsSide(glob_rel))) {    switch(RelationTermType(Infix_RhsSide(glob_rel))) {
3362    case e_real:    case e_real:
3363    case e_int:    case e_int:
# Line 3350  static void SetUpInvertTokenTop(struct r Line 3381  static void SetUpInvertTokenTop(struct r
3381  }  }
3382    
3383  /**  /**
3384   *** RelationInvertToken inverts tokens until the variable      @see RelationFindRoots
  *** being solved for is found.  It is assumed that this  
  *** variable only resides at ONE leaf of the relation tree.  
  *** It is the calling function's responsibility to make sure  
  *** this is the case and call another solver if needed.  
  *** If the variable (with varnum = glob_varnum) is found,  
  *** the solution list will contain all solutions to the  
  *** equations.  It is the calling function's responsibility  
  *** to select the root that suits his needs.  
  *** Returns TRUE for success, FALSE for failure.  
  **/  
3385    
3386  /* Note that there appears to be some redundant checking here and      Invert tokens until the variable being solved for is found.  
3387   * we could probably be more efficient  
3388   */      It is assumed that this
3389        variable only resides at ONE leaf of the relation tree.
3390        It is the calling function's responsibility to make sure
3391        this is the case and call another solver if needed.
3392        If the variable (with varnum = glob_varnum) is found,
3393        the solution list will contain all solutions to the
3394        equations.  It is the calling function's responsibility
3395        to select the root that suits his needs.
3396        Returns TRUE for success, FALSE for failure.
3397    
3398        @NOTE Note that there appears to be some redundant checking here and
3399        we could probably be more efficient
3400    */
3401  int RelationInvertToken(struct relation_term **term,  int RelationInvertToken(struct relation_term **term,
3402                  struct ds_soln_list *soln_list,                  struct ds_soln_list *soln_list,
3403                  enum safe_err *not_safe)                  enum safe_err *not_safe
3404  {  ){
3405    int side,ndx;    int side,ndx;
3406    double value = 0.0;    double value = 0.0;
3407    struct relation_term *invert_side;    struct relation_term *invert_side;
# Line 3665  int RelationInvertToken(struct relation_ Line 3698  int RelationInvertToken(struct relation_
3698  }  }
3699    
3700  /*  /*
3701   * RelationInvertTokenTop is baisically a while loop which      @see RelationFindRoots
3702   * calls RelationInvertToken.  See RelationInvertToken for  
3703   * information on what this function does.      RelationInvertTokenTop is baisically a while loop which
3704   */      calls RelationInvertToken.  See RelationInvertToken for
3705  int RelationInvertTokenTop(struct ds_soln_list *soln_list)      information on what this function does.
3706  {  */
3707    int RelationInvertTokenTop(struct ds_soln_list *soln_list){
3708    int result;    int result;
3709    struct relation_term *invert_side;    struct relation_term *invert_side;
3710    enum safe_err not_safe = safe_ok;    enum safe_err not_safe = safe_ok;
# Line 3687  int RelationInvertTokenTop(struct ds_sol Line 3721  int RelationInvertTokenTop(struct ds_sol
3721    return result;    return result;
3722  }  }
3723    
3724  /*************************************************************************/  /*-----------------------------------------------------------------------------
3725  /*************************Rootfinding Functions***************************/    ROOT-FINDING FUNCTIONS
3726  /*************************************************************************/    to support the RelationFindRoots function.
3727    */
3728    
3729  /*  /**
3730   *  --- glob_rel is ASSUMED to be of type e_token. ---      Set the value of the variable being solved
3731   *  CalcResidGivenValue sets the value of the variable being solved      for (given the varnum) and calculate the residual.  
3732   *  for (given the varnum) and calculates the residual.  Note that  
3733   *  this function uses the glob_rel which should have been set in      @NOTE glob_rel is ASSUMED to be of type e_token. ---
3734   *  RelationFindRoots (and reduced by SearchEval_Branch).  This  
3735   *  functions takes an excessive number of arguments so it will      @NOTE uses the glob_rel which should have been set in
3736   *  look like an ExtEvalFunc to our rootfinder.      RelationFindRoots (and reduced by SearchEval_Branch).  This
3737   */      functions takes an excessive number of arguments so it will
3738        look like an ExtEvalFunc to our root-finder.
3739    */
3740  static  static
3741  int CalcResidGivenValue(int *mode, int *m, unsigned long *varnum,  int CalcResidGivenValue(int *mode, int *m, unsigned long *varnum,
3742                  double *val, double *u, double *f, double *g)                  double *val, double *u, double *f, double *g)
# Line 3710  int CalcResidGivenValue(int *mode, int * Line 3747  int CalcResidGivenValue(int *mode, int *
3747     *  glob_rel is ASSUMED to be of type e_token.     *  glob_rel is ASSUMED to be of type e_token.
3748     */     */
3749    
3750    (void)mode;  /*  stop gcc whine about unused parameter  */    UNUSED_PARAMETER(mode);
3751    (void)u;     /*  stop gcc whine about unused parameter  */    UNUSED_PARAMETER(u);
3752    (void)g;     /*  stop gcc whine about unused parameter  */    UNUSED_PARAMETER(g);
3753    
3754    SetRealAtomValue(    SetRealAtomValue(
3755        ((struct Instance *)gl_fetch(RelationVarList(glob_rel),*varnum)),        ((struct Instance *)gl_fetch(RelationVarList(glob_rel),*varnum)),
# Line 3739  int CalcResidGivenValue(int *mode, int * Line 3776  int CalcResidGivenValue(int *mode, int *
3776    return 0;    return 0;
3777  }  }
3778    
3779  /*  /**
3780   *  RootFind is a distributor to a rootfinding method zbrent.      RootFind is a distributor to a root-finding method zbrent.
3781   *  at present, the rootfind can only handle token relations.      at present, the root-find can only handle token relations.
3782   * if varnum is 0 and status is NULL, frees the internal      if varnum is 0 and status is NULL, frees the internal
3783   * memory recycle and returns. Currently this is slaved      memory recycle and returns. Currently this is slaved
3784   * to RelationFindRoots internal reset.      to RelationFindRoots internal reset.
3785   * This function is not threadsafe.      This function is not threadsafe.
3786   */  */
3787  static  static
3788  double RootFind(struct relation *rel,  double RootFind(struct relation *rel,
3789              double *lower_bound,          double *lower_bound, double *upper_bound,
3790              double *upper_bound,          double *nominal,
3791              double *nominal,          double *tolerance,
3792              double *tolerance,          unsigned long varnum,
3793              unsigned long varnum,          int *status
3794                  int *status)  ){
 {  
3795    double *f = NULL; /* vector of residuals, borrowed from tmpalloc */    double *f = NULL; /* vector of residuals, borrowed from tmpalloc */
3796    ExtEvalFunc *func;    ExtEvalFunc *func;
3797    int mode;             /* to pass to the eval func */    int mode;             /* to pass to the eval func */
# Line 3801  double RootFind(struct relation *rel, Line 3837  double RootFind(struct relation *rel,
3837                  &(m),&(n),x,u,f,g,tolerance,status);                  &(m),&(n),x,u,f,g,tolerance,status);
3838  }  }
3839    
3840  /*************************************************************************/  /*------------------------------------------------------------------------------
3841  /*************************External Function******************************/    UTILITY FUNCTIONS
3842  /*************************************************************************/  */
   
 /**  
  *** RelationFindRoot WILL find a root if there is one. It is  
  *** in charge of trying every trick in the book. It returns  
  *** 1 for success and 0 for failure. In general compiler functions  
  *** return 0 for success but this function returns 1 for success  
  *** because success = 1 is the convention on the solver side.  
  *** (we really should make a system wide convention)  
  *** A return of -1 indicates a problem such as var not found.  
  *** If nsolns > 1 then a list of solutions will be returned.  
  **/  
3843    
3844  /* Note we should recycle the memory used for glob_rel */  /**
3845  double *RelationFindRoots(struct Instance *i,      Create a static buffer to use for temporary memory storage
                   double lower_bound,  
                   double upper_bound,  
                   double nominal,  
                   double tolerance,  
                   unsigned long *varnum,  
                   int *able,  
                   int *nsolns)  
 {  
   struct relation *rel;  
   double sideval;  
   enum Expr_enum reltype;  
   static struct ds_soln_list soln_list = {0,0,NULL};  
   CONST struct gl_list_t *list;  
3846    
3847    /* check for recycle shutdown */      Temporarily allocates a given number of bytes.  The memory need
3848    if (i==NULL && varnum == NULL && able == NULL && nsolns == NULL) {      not be freed, but the next call to this function will reuse the
3849      if (soln_list.soln != NULL) {      previous allocation. Memory returned will NOT be zeroed.
3850        ascfree(soln_list.soln);      Calling with nbytes==0 will free any memory allocated.
3851        soln_list.soln = NULL;  */
3852        soln_list.length = soln_list.capacity = 0;  char *tmpalloc(int nbytes){
3853      static char *ptr = NULL;
3854      static int cap = 0;
3855    
3856      if( nbytes > 0 ) {
3857        if( nbytes > cap ) {
3858          if( ptr ) ascfree(ptr);
3859          ptr = (char *)ascmalloc(nbytes);
3860          cap = nbytes;
3861      }      }
     RootFind(NULL,NULL,NULL,NULL,NULL,0L,NULL); /*clear brent recycle */  
     RelationCreateTmp(0,0,e_nop); /* clear tmprelation recycle */  
     return NULL;  
3862    }    }
3863    /* check assertions */    else {
3864  #ifndef NDEBUG      if( ptr ) ascfree(ptr);
3865    if( i == NULL ) {      ptr = NULL;
3866      FPRINTF(ASCERR, "error in RelationFindRoot: NULL instance\n");      cap = 0;
     glob_rel = NULL;  
     return NULL;  
   }  
   if (able == NULL){  
     FPRINTF(ASCERR,"error in RelationFindRoot: NULL able ptr\n");  
     glob_rel = NULL;  
     return NULL;  
   }  
   if (varnum == NULL){  
     FPRINTF(ASCERR,"error in RelationFindRoot: NULL varnum\n");  
     glob_rel = NULL;  
     return NULL;  
   }  
   if( InstanceKind(i) != REL_INST ) {  
     FPRINTF(ASCERR, "error in RelationFindRoot: not relation\n");  
     glob_rel = NULL;  
     return NULL;  
3867    }    }
3868  #endif    return ptr;
3869    }
3870    
3871    *able = FALSE;  /**
3872    *nsolns = -1;     /* nsolns will be -1 for a very unhappy rootfinder */      @TODO what this does needs to be documented here
   glob_rel = NULL;  
   glob_done = 0;  
   soln_list.length = 0; /* reset len to 0. if NULL to start, append mallocs */  
   append_soln(&soln_list,0.0);  
   rel = (struct relation *)GetInstanceRelation(i, &reltype);  
   if( rel == NULL ) {  
     FPRINTF(ASCERR, "error in RelationFindRoot: NULL relation\n");  
     glob_rel = NULL; return NULL;  
   }  
   /* here we should switch and handle all types. at present we don't  
    * handle anything except e_token  
    */  
   if( reltype != e_token ) {  
     FPRINTF(ASCERR, "error in RelationFindRoot: non-token relation\n");  
     glob_rel = NULL;  
     return NULL;  
   }  
3873    
3874    if (RelationRelop(rel) == e_equal){      it is questionable whether this should be unified with the
3875      glob_rel = RelationTmpTokenCopy(rel);      ArgsForToken function in relation.c
3876      assert(glob_rel!=NULL);  */
3877      glob_done = 0;  int ArgsForRealToken(enum Expr_enum type){
3878      list = RelationVarList(glob_rel);     switch(type) {
3879      if( *varnum >= 1 && *varnum <= gl_length(list)){     case e_zero:
3880        glob_done = 1;     case e_int:
3881      }     case e_real:
3882      if (!glob_done) {     case e_var:
3883        FPRINTF(ASCERR, "error in FindRoot: var not found\n");       return(0);
       glob_rel = NULL;  
       return NULL;  
     }  
3884    
3885      glob_varnum = *varnum;     case e_func:
3886      glob_done = 0;     case e_uminus:
3887      assert(Infix_LhsSide(glob_rel) != NULL);        return(1);
     /* In the following if statements we look for the target variable  
      * to the left and right, evaluating all branches without the  
      * target.  
      */  
     if (SearchEval_Branch(Infix_LhsSide(glob_rel)) < 1) {  
       sideval = RelationBranchEvaluator(Infix_LhsSide(glob_rel));  
       if (finite(sideval)) {  
         InsertBranchResult(Infix_LhsSide(glob_rel),sideval);  
       } else {  
         FPRINTF(ASCERR,"Inequality in RelationFindRoots. Infinite RHS.\n");  
         glob_rel = NULL;  
         return NULL;  
       }  
     }  
     assert(Infix_RhsSide(glob_rel) != NULL);  
     if (SearchEval_Branch(Infix_RhsSide(glob_rel)) < 1) {  
         sideval = RelationBranchEvaluator(Infix_RhsSide(glob_rel));  
         if (finite(sideval)) {  
           InsertBranchResult(Infix_RhsSide(glob_rel),sideval);  
         } else {  
           FPRINTF(ASCERR,"Inequality in RelationFindRoots. Infinite LHS.\n");  
           glob_rel = NULL;  
           return NULL;  
         }  
     }  
     if (glob_done < 1) {  
       /* RelationInvertToken never found variable */  
       glob_done = 0;  
       *able = FALSE;  
       return soln_list.soln;  
     }  
     if (glob_done == 1) {  
       /* set to 0 so while loop in RelationInvertToken will work */  
       glob_done = 0;  
       glob_done = RelationInvertTokenTop(&(soln_list));  
     }  
     if (glob_done == 1) { /* if still one, token inversions successful */  
       glob_done = 0;  
       *nsolns= soln_list.length;  
       *able = TRUE;  
       return soln_list.soln;  
     }  
     /* CALL ITERATIVE SOLVER */  
     *soln_list.soln = RootFind(glob_rel,&(lower_bound),  
                        &(upper_bound),&(nominal),  
                        &(tolerance),  
                        glob_varnum,able);  
3888    
3889      glob_done = 0;     case e_plus:
3890      if(*able == 0) { /* Root-Find returns 0 for success*/     case e_minus:
3891        *nsolns = 1;     case e_times:
3892        *able = TRUE;     case e_divide:
3893      } else {     case e_power:
3894        *able = FALSE;     case e_ipower:
3895      }        return(2);
     return soln_list.soln;  
3896    
3897    }     default:
3898    FPRINTF(ASCERR,"Inequality in RelationFindRoots. can't find roots.\n");        FPRINTF(ASCERR,"ArgsForRealToken: Unknown relation term type.\n");
3899    *able = FALSE;        return(0);
3900    return soln_list.soln;     }
3901  }  }
3902    
3903    static int IsZero(struct dimnode *node){
3904        if( node->type==e_zero || (node->type==e_real && node->real_const==0.0) ){
3905            return TRUE;
3906        }
3907        return FALSE;
3908    }
3909    
3910    
3911    #define START 10000  /* largest power of 10 held by a short */
3912    
3913  /*  static struct fraction real_to_frac(double real){
3914   *  Temporary Functions for testing direct solve.     short  num, den;
3915   *  Remove calls from interface.c when this is removed.     for( den=START; den>1 && fabs(real)*den>SHRT_MAX; den /= 10 ) ;
3916   */     num = (short)(fabs(real)*den + 0.5);
3917  void PrintDirectResult(struct Instance *i)     if( real < 0.0 ) num = -num;
3918  {     return( CreateFraction(num,den) );
3919    }
3920    
3921    #undef START
3922    
3923    
3924    /*------------------------------------------------------------------------------
3925        Temporary Functions for testing direct solve.
3926        Remove calls from interface.c when this is removed.
3927    */
3928    
3929    void PrintDirectResult(struct Instance *i){
3930    struct relation *rel;    struct relation *rel;
3931    enum Expr_enum reltype;    enum Expr_enum reltype;
3932    int num,status,n,nsoln;    int num,status,n,nsoln;
# Line 4037  void CollectShares(struct Instance *i,st Line 3988  void CollectShares(struct Instance *i,st
3988    }    }
3989  }  }
3990    
3991    
3992  struct gl_list_t *  struct gl_list_t *
3993  CollectTokenRelationsWithUniqueBINlessShares(struct Instance *i,  CollectTokenRelationsWithUniqueBINlessShares(struct Instance *i,
3994                                               unsigned long maxlen)                                               unsigned long maxlen)
# Line 4063  CollectTokenRelationsWithUniqueBINlessSh Line 4015  CollectTokenRelationsWithUniqueBINlessSh
4015      return data.list;      return data.list;
4016    }    }
4017  }  }
4018    
4019    #ifndef NDEBUG
4020    /**
4021        Utility function to perform debug checking of (input) instance and residual
4022        (or gradient) (output) pointers in the various functions in this file.
4023    */
4024    static int  relutil_check_inst_and_res(struct Instance *i, double *res){
4025    # ifdef RELUTIL_CHECK_ABORT
4026        if(i==NULL){
4027            Asc_Panic(2,__FUNCTION__,"NULL instance");
4028        }else if (res==NULL){
4029            Asc_Panic(2,__FUNCTION__,"NULL residual pointer");
4030        }else if(InstanceKind(i)!=REL_INST){
4031            Asc_Panic(2,__FUNCTION__,"Not a relation");
4032        }
4033    # else
4034      if( i == NULL ) {
4035        ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL instance");
4036        return 0;
4037      }else if (res == NULL){
4038        ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL residual ptr");
4039        return 0;
4040      }else if( InstanceKind(i) != REL_INST ) {
4041        ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation");
4042        return 0;
4043      }
4044      return 1;
4045    # endif
4046    }
4047    #endif  

Legend:
Removed from v.693  
changed lines
  Added in v.694

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