/[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 1128 by johnpye, Sat Jan 13 08:33:43 2007 UTC revision 1129 by johnpye, Sat Jan 13 11:40:59 2007 UTC
# Line 1  Line 1 
1  /*  ASCEND modelling environment  /*  ASCEND modelling environment
2      Copyright 1997, Carnegie Mellon University      Copyright 1997, Carnegie Mellon University
3      Copyright (C) 2006 Carnegie Mellon University      Copyright (C) 2006-2007 Carnegie Mellon University
4    
5      This program is free software; you can redistribute it and/or modify      This program is free software; you can redistribute it and/or modify
6      it under the terms of the GNU General Public License as published by      it under the terms of the GNU General Public License as published by
# Line 16  Line 16 
16      along with this program; if not, write to the Free Software      along with this program; if not, write to the Free Software
17      Foundation, Inc., 59 Temple Place - Suite 330,      Foundation, Inc., 59 Temple Place - Suite 330,
18      Boston, MA 02111-1307, USA.      Boston, MA 02111-1307, USA.
19  *//**  *//** @file
20      @file      LSODE Integrator
21      LSODE integrator.  
22        Here we are solving the system of first order ordinary differential
23        equations, ydot = F(y,t), with y(t_0)=y_0 given.
24    
25        In some places 'x' is named instead of 't'. We're trying to migrate to using
26        't' to label the independent variable.
27    
28        There is some excellent documentation about the LSODE integrator available
29        in the document "Description and Use of LSODE, the Livermore Solver for
30        Ordinary Differential Equations, by K Radhakrishnan and A C Hindmarsh.
31        http://www.llnl.gov/CASC/nsde/pubs/u113855.pdf
32    
33      (old implementation notes:)      (old implementation notes:)
34    
# Line 34  Line 44 
44      a solve.      a solve.
45    
46      @NOTE The above doesn't work since lsode doesn't use the same t internally      @NOTE The above doesn't work since lsode doesn't use the same t internally
47      that we hand it.      that we hand it. @ENDNOTE
   
48  *//*  *//*
49      by Kirk Abbott and Ben Allan      original version by Kirk Abbott and Ben Allan, Jan 1994.
50      Created: 1/94      Last in CVS $Revision: 1.29 $ $Date: 2000/01/25 02:26:31 $ $Author: ballan $
51      Version: $Revision: 1.29 $  */
     Version control file: $RCSfile: Lsode.c,v $  
     Date last modified: $Date: 2000/01/25 02:26:31 $  
     Last modified by: $Author: ballan $  
  */  
52    
53  #ifndef NO_SIGNAL_TRAPS  #ifndef NO_SIGNAL_TRAPS
54  #include <signal.h>  #include <signal.h>
# Line 81  const IntegratorInternals integrator_lso Line 86  const IntegratorInternals integrator_lso
86      ,integrator_lsode_params_default      ,integrator_lsode_params_default
87      ,integrator_analyse_ode /* note, this routine is back in integrator.c */      ,integrator_analyse_ode /* note, this routine is back in integrator.c */
88      ,integrator_lsode_solve      ,integrator_lsode_solve
89        ,integrator_lsode_write_matrix
90      ,integrator_lsode_free      ,integrator_lsode_free
91      ,INTEG_LSODE      ,INTEG_LSODE
92      ,"LSODE"      ,"LSODE"
# Line 182  typedef struct IntegratorLsodeDataStruct Line 188  typedef struct IntegratorLsodeDataStruct
188      long n_eqns;                     /**< dimension of state vector */      long n_eqns;                     /**< dimension of state vector */
189      int *input_indices;              /**< vector of state vars indexes */      int *input_indices;              /**< vector of state vars indexes */
190      int *output_indices;             /**< vector of derivative var indexes */      int *output_indices;             /**< vector of derivative var indexes */
191      struct var_variable **y_vars;    /**< NULL terminated list of states vars */      struct var_variable **y_vars;    /**< NULL-terminated list of states vars */
192      struct var_variable **ydot_vars; /**< NULL terminated list of derivative vars*/      struct var_variable **ydot_vars; /**< NULL-terminated list of derivative vars*/
193      struct rel_relation **rlist;     /**< NULL terminated list of relevant rels      struct rel_relation **rlist;     /**< NULL-terminated list of relevant rels
194                                          to be differentiated */                                          to be differentiated */
195      DenseMatrix dydx_dx;             /**< change in derivatives wrt states      DenseMatrix dydot_dy;               /**< change in derivatives wrt states */
                                         I prefer to call this: d(ydot)/dy */  
196    
197      IntegratorLsodeLastCallType lastcall;  /* type of last call; func or grad */      IntegratorLsodeLastCallType lastcall;  /* type of last call; func or grad */
198      IntegratorLsodeStatusCode   status;    /* solve status */      IntegratorLsodeStatusCode   status;    /* solve status */
# Line 269  void integrator_lsode_create(IntegratorS Line 274  void integrator_lsode_create(IntegratorS
274      d->y_vars=NULL;      d->y_vars=NULL;
275      d->ydot_vars=NULL;      d->ydot_vars=NULL;
276      d->rlist=NULL;      d->rlist=NULL;
277      d->dydx_dx=(DenseMatrix){NULL,0,0};      d->dydot_dy=DENSEMATRIX_EMPTY;
278      blsys->enginedata=(void*)d;      blsys->enginedata=(void*)d;
279      integrator_lsode_params_default(blsys);      integrator_lsode_params_default(blsys);
280    
# Line 297  void integrator_lsode_free(void *engined Line 302  void integrator_lsode_free(void *engined
302      if(d.rlist)ASC_FREE(d.rlist);      if(d.rlist)ASC_FREE(d.rlist);
303      d.rlist =  NULL;      d.rlist =  NULL;
304    
305      densematrix_destroy(d.dydx_dx);      densematrix_destroy(d.dydot_dy);
306    
307      d.n_eqns = 0L;      d.n_eqns = 0L;
308  }  }
# Line 620  int integrator_lsode_derivatives(Integra Line 625  int integrator_lsode_derivatives(Integra
625    asc_assert(blsys!=NULL);    asc_assert(blsys!=NULL);
626    enginedata = (IntegratorLsodeData *)blsys->enginedata;    enginedata = (IntegratorLsodeData *)blsys->enginedata;
627    asc_assert(enginedata!=NULL);    asc_assert(enginedata!=NULL);
628    asc_assert(DENSEMATRIX_DATA(enginedata->dydx_dx)!=NULL);    asc_assert(DENSEMATRIX_DATA(enginedata->dydot_dy)!=NULL);
629    asc_assert(enginedata->input_indices!=NULL);    asc_assert(enginedata->input_indices!=NULL);
630    
631    double **dy_dx = DENSEMATRIX_DATA(enginedata->dydx_dx);    double **dydot_dy = DENSEMATRIX_DATA(enginedata->dydot_dy);
632    int *inputs_ndx_list = enginedata->input_indices;    int *inputs_ndx_list = enginedata->input_indices;
633    int *outputs_ndx_list = enginedata->output_indices;    int *outputs_ndx_list = enginedata->output_indices;
634    asc_assert(ninputs == blsys->n_y);    asc_assert(ninputs == blsys->n_y);
# Line 660  int integrator_lsode_derivatives(Integra Line 665  int integrator_lsode_derivatives(Integra
665      FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");      FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
666      goto error;      goto error;
667    }    }
668    result = Compute_dy_dx_smart(blsys->system, scratch_vector, dy_dx,    result = Compute_dy_dx_smart(blsys->system, scratch_vector, dydot_dy,
669                                 inputs_ndx_list, ninputs,                                 inputs_ndx_list, ninputs,
670                                 outputs_ndx_list, noutputs);                                 outputs_ndx_list, noutputs);
671    
# Line 817  static void LSODE_JEX(int *neq ,double * Line 822  static void LSODE_JEX(int *neq ,double *
822      Map data from C based matrix to Fortan matrix.      Map data from C based matrix to Fortan matrix.
823      We will send in a column major ordering vector for pd.      We will send in a column major ordering vector for pd.
824    */    */
825    asc_assert(*neq == DENSEMATRIX_NCOLS(lsodedata->dydx_dx));    asc_assert(*neq == DENSEMATRIX_NCOLS(lsodedata->dydot_dy));
826    asc_assert(*nrpd == DENSEMATRIX_NROWS(lsodedata->dydx_dx));    asc_assert(*nrpd == DENSEMATRIX_NROWS(lsodedata->dydot_dy));
827    for (j=0;j<*neq;j++) { /* loop through columnns */    for (j=0;j<*neq;j++) { /* loop through columnns */
828      for (i=0;i<*nrpd;i++){ /* loop through rows */      for (i=0;i<*nrpd;i++){ /* loop through rows */
829        /* 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.dydot_dy[i][j]); */
830        *pd++ = DENSEMATRIX_ELEM(lsodedata->dydx_dx,i,j);        *pd++ = DENSEMATRIX_ELEM(lsodedata->dydot_dy,i,j);
831      }      }
832    }    }
833    
# Line 867  int integrator_lsode_solve(IntegratorSys Line 872  int integrator_lsode_solve(IntegratorSys
872    
873      d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);      d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
874      d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);      d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
875      d->dydx_dx = densematrix_create(d->n_eqns,d->n_eqns);      d->dydot_dy = densematrix_create(d->n_eqns,d->n_eqns);
876    
877      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);
878      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);
# Line 1207  void XASCWV( char *msg, /* pointer to st Line 1212  void XASCWV( char *msg, /* pointer to st
1212      }      }
1213      error_reporter_end_flush();      error_reporter_end_flush();
1214  }  }
1215    
1216    int integrator_lsode_write_matrix(const IntegratorSystem *blsys, FILE *fp){
1217        IntegratorLsodeData *enginedata;
1218        asc_assert(blsys!=NULL);
1219        asc_assert(blsys->engine==INTEG_LSODE);
1220        asc_assert(blsys->enginedata);
1221        enginedata = (IntegratorLsodeData *)blsys->enginedata;
1222    
1223        if(!DENSEMATRIX_DATA(enginedata->dydot_dy)){
1224            ERROR_REPORTER_HERE(ASC_PROG_ERR,"dydot_dy contains no data");
1225        }
1226    
1227        densematrix_write_mmio(enginedata->dydot_dy,fp);
1228        CONSOLE_DEBUG("Returning after matrix output");
1229        return 0;
1230    }
1231        
1232    
1233    
1234    
1235    
Legend:
Removed from v.1128  
changed lines
  Added in v.1129
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22