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

Diff of /trunk/base/generic/solver/lsode.c

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

revision 1127 by johnpye, Wed Jan 10 06:26:02 2007 UTC revision 1128 by johnpye, Sat Jan 13 08:33:43 2007 UTC
# Line 71  Line 71 
71  #include "slv_client.h"  #include "slv_client.h"
72  #include "slv_stdcalls.h"  #include "slv_stdcalls.h"
73  #include <packages/sensitivity.h>  #include <packages/sensitivity.h>
74    #include "densemtx.h"
75    
76  #include "integrator.h"  #include "integrator.h"
77  #include "lsode.h"  #include "lsode.h"
# Line 177  typedef enum{ Line 178  typedef enum{
178      Enumeration to tell ASCEND if anything failed in a FEX or JEX call.      Enumeration to tell ASCEND if anything failed in a FEX or JEX call.
179  */  */
180    
181  typedef struct{  typedef struct IntegratorLsodeDataStruct{
     IntegratorLsodeLastCallType lastcall;  /* type of last call; func or grad */  
     IntegratorLsodeStatusCode   status;    /* solve status */  
     char stop;                             /* stop requested? */  
     int partitioned;                       /* partioned func evals or not */  
 } IntegratorLsodeStatus;  
 /**<  
     Bits of data that LSODE needs to be able to send/recieve from inside a  
     JEX or FEX call.  
 */  
   
     typedef struct IntegratorLsodeDataStruct{  
182      long n_eqns;                     /**< dimension of state vector */      long n_eqns;                     /**< dimension of state vector */
183      int *input_indices;              /**< vector of state vars indexes */      int *input_indices;              /**< vector of state vars indexes */
184      int *output_indices;             /**< vector of derivative var indexes */      int *output_indices;             /**< vector of derivative var indexes */
# Line 196  typedef struct{ Line 186  typedef struct{
186      struct var_variable **ydot_vars; /**< NULL terminated list of derivative vars*/      struct var_variable **ydot_vars; /**< NULL terminated list of derivative vars*/
187      struct rel_relation **rlist;     /**< NULL terminated list of relevant rels      struct rel_relation **rlist;     /**< NULL terminated list of relevant rels
188                                          to be differentiated */                                          to be differentiated */
189      double **dydx_dx;                /**< change in derivatives wrt states      DenseMatrix dydx_dx;             /**< change in derivatives wrt states
190                                          I prefer to call this: d(ydot)/dy */                                          I prefer to call this: d(ydot)/dy */
191    
192      IntegratorLsodeLastCallType lastcall;  /* type of last call; func or grad */      IntegratorLsodeLastCallType lastcall;  /* type of last call; func or grad */
# Line 248  typedef void LsodeJacobianFn(int *, doub Line 238  typedef void LsodeJacobianFn(int *, doub
238  */  */
239    
240  int integrator_lsode_setup_diffs(IntegratorSystem *blsys);  int integrator_lsode_setup_diffs(IntegratorSystem *blsys);
 static double **lsode_densematrix_create(int nrows, int ncols);  
 static void lsode_densematrix_destroy(double **matrix,int nrows);  
241    
242  /**  /**
243      void LSODE(&fex, &neq, y, &x, &xend, &itol, reltol, abtol, &itask,      void LSODE(&fex, &neq, y, &x, &xend, &itol, reltol, abtol, &itask,
# Line 281  void integrator_lsode_create(IntegratorS Line 269  void integrator_lsode_create(IntegratorS
269      d->y_vars=NULL;      d->y_vars=NULL;
270      d->ydot_vars=NULL;      d->ydot_vars=NULL;
271      d->rlist=NULL;      d->rlist=NULL;
272      d->dydx_dx=NULL;      d->dydx_dx=(DenseMatrix){NULL,0,0};
273      blsys->enginedata=(void*)d;      blsys->enginedata=(void*)d;
274      integrator_lsode_params_default(blsys);      integrator_lsode_params_default(blsys);
275    
# Line 309  void integrator_lsode_free(void *engined Line 297  void integrator_lsode_free(void *engined
297      if(d.rlist)ASC_FREE(d.rlist);      if(d.rlist)ASC_FREE(d.rlist);
298      d.rlist =  NULL;      d.rlist =  NULL;
299    
300      if(d.dydx_dx!=NULL){      densematrix_destroy(d.dydx_dx);
         lsode_densematrix_destroy(d.dydx_dx, d.n_eqns);  
         d.dydx_dx =  NULL;  
     }  
301    
302      d.n_eqns = 0L;      d.n_eqns = 0L;
303  }  }
# Line 402  int integrator_lsode_params_default(Inte Line 387  int integrator_lsode_params_default(Inte
387      return 0;      return 0;
388  }    }  
389    
 /*---------------------------------------------------------  
   Couple of matrix methods...?  
 */  
   
 static double **lsode_densematrix_create(int nrows, int ncols){  
   int c;  
   double **result;  
   assert(nrows>0);  
   assert(ncols>0);  
   result = ASC_NEW_ARRAY(double *, nrows);  
   for (c=0;c<nrows;c++) {  
     result[c] = ASC_NEW_ARRAY_CLEAR(double, ncols);  
   }  
   return result;  
 }  
   
 static void lsode_densematrix_destroy(double **matrix,int nrows){  
   int c;  
   if (matrix) {  
     for (c=0;c<nrows;c++) {  
       if (matrix[c]) {  
         ascfree((char *)matrix[c]);  
       }  
     }  
     ascfree((char *)matrix);  
   }  
 }  
   
390  /*------------------------------------------------------------------------------  /*------------------------------------------------------------------------------
391    PROBLEM ANALYSIS    PROBLEM ANALYSIS
392  */  */
# Line 663  int integrator_lsode_derivatives(Integra Line 620  int integrator_lsode_derivatives(Integra
620    asc_assert(blsys!=NULL);    asc_assert(blsys!=NULL);
621    enginedata = (IntegratorLsodeData *)blsys->enginedata;    enginedata = (IntegratorLsodeData *)blsys->enginedata;
622    asc_assert(enginedata!=NULL);    asc_assert(enginedata!=NULL);
623    asc_assert(enginedata->dydx_dx!=NULL);    asc_assert(DENSEMATRIX_DATA(enginedata->dydx_dx)!=NULL);
624    asc_assert(enginedata->input_indices!=NULL);    asc_assert(enginedata->input_indices!=NULL);
625    
626    double **dy_dx = enginedata->dydx_dx;    double **dy_dx = DENSEMATRIX_DATA(enginedata->dydx_dx);
627    int *inputs_ndx_list = enginedata->input_indices;    int *inputs_ndx_list = enginedata->input_indices;
628    int *outputs_ndx_list = enginedata->output_indices;    int *outputs_ndx_list = enginedata->output_indices;
629    asc_assert(ninputs == blsys->n_y);    asc_assert(ninputs == blsys->n_y);
# Line 860  static void LSODE_JEX(int *neq ,double * Line 817  static void LSODE_JEX(int *neq ,double *
817      Map data from C based matrix to Fortan matrix.      Map data from C based matrix to Fortan matrix.
818      We will send in a column major ordering vector for pd.      We will send in a column major ordering vector for pd.
819    */    */
820      asc_assert(*neq == DENSEMATRIX_NCOLS(lsodedata->dydx_dx));
821      asc_assert(*nrpd == DENSEMATRIX_NROWS(lsodedata->dydx_dx));
822    for (j=0;j<*neq;j++) { /* loop through columnns */    for (j=0;j<*neq;j++) { /* loop through columnns */
823      for (i=0;i<*nrpd;i++){ /* loop through rows */      for (i=0;i<*nrpd;i++){ /* loop through rows */
824        /* CONSOLE_DEBUG("JAC[r=%d,c=%d]=%f",i,j,lsodedata.dydx_dx[i][j]); */        /* CONSOLE_DEBUG("JAC[r=%d,c=%d]=%f",i,j,lsodedata.dydx_dx[i][j]); */
825        *pd++ = lsodedata->dydx_dx[i][j];        *pd++ = DENSEMATRIX_ELEM(lsodedata->dydx_dx,i,j);
826      }      }
827    }    }
828    
# Line 908  int integrator_lsode_solve(IntegratorSys Line 867  int integrator_lsode_solve(IntegratorSys
867    
868      d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);      d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
869      d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);      d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
870      d->dydx_dx = lsode_densematrix_create(d->n_eqns,d->n_eqns);      d->dydx_dx = densematrix_create(d->n_eqns,d->n_eqns);
871    
872      d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);      d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);
873      d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);      d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);
Legend:
Removed from v.1127  
changed lines
  Added in v.1128
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22