/[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 965 by johnpye, Tue Dec 12 13:36:52 2006 UTC revision 966 by johnpye, Thu Dec 14 14:04:54 2006 UTC
# Line 814  static void LSODE_JEX(int *neq ,double * Line 814  static void LSODE_JEX(int *neq ,double *
814    
815  /**  /**
816      The public function: here we do the actual integration, I guess.      The public function: here we do the actual integration, I guess.
817    
818        Return 1 on success
819  */  */
820  int integrator_lsode_solve(IntegratorSystem *blsys  int integrator_lsode_solve(IntegratorSystem *blsys
821          , unsigned long start_index, unsigned long finish_index          , unsigned long start_index, unsigned long finish_index
# Line 886  int integrator_lsode_solve(IntegratorSys Line 888  int integrator_lsode_solve(IntegratorSys
888    neq = blsys->n_y;    neq = blsys->n_y;
889    nobs = blsys->n_obs;    nobs = blsys->n_obs;
890    
891    x[0] = integrator_get_t(blsys);    /* samplelist_debug(blsys->samples); */
892    
893      /* x[0] = integrator_get_t(blsys); */
894      x[0] = integrator_getsample(blsys, 0);
895    x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */    x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */
896    lrw = 22 + 9*neq + neq*neq;    lrw = 22 + 9*neq + neq*neq;
897    rwork = ASC_NEW_ARRAY_CLEAR(double, lrw+1);    rwork = ASC_NEW_ARRAY_CLEAR(double, lrw+1);
# Line 917  int integrator_lsode_solve(IntegratorSys Line 922  int integrator_lsode_solve(IntegratorSys
922    iwork[5] = integrator_get_maxsubsteps(blsys);    iwork[5] = integrator_get_maxsubsteps(blsys);
923    mf = 21;      /* 21 = BDF with exact jacobian. 22 = BDF with finite diff Jacobian */    mf = 21;      /* 21 = BDF with exact jacobian. 22 = BDF with finite diff Jacobian */
924    
925      if(x[0] > integrator_getsample(blsys, 2)){
926        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value");
927        return 0;
928      }
929    
930    /* put the values from derivative system into the record */    /* put the values from derivative system into the record */
931    integrator_setsample(blsys, start_index, x[0]);    integrator_setsample(blsys, start_index, x[0]);
932    
# Line 930  int integrator_lsode_solve(IntegratorSys Line 940  int integrator_lsode_solve(IntegratorSys
940      the loop ahead in time, all we need to do is keep upping      the loop ahead in time, all we need to do is keep upping
941      xend.      xend.
942    */    */
943      
944    blsys->currentstep = 0;    blsys->currentstep = 0;
945    for (index = start_index; index < finish_index; index++,  blsys->currentstep++) {    for (index = start_index; index < finish_index; index++,  blsys->currentstep++) {
946      xend = integrator_getsample(blsys, index+1);      xend = integrator_getsample(blsys, index+1);
947      xprev = x[0];      xprev = x[0];
948      /* CONSOLE_DEBUG("BEFORE %lu LSODE CALL\n", index); */      asc_assert(xend > xprev);
949        /* CONSOLE_DEBUG("LSODE call #%lu: x = [%f,%f]", index,xprev,xend); */
950    
951  # ifndef NO_SIGNAL_TRAPS  # ifndef NO_SIGNAL_TRAPS
952      if (setjmp(g_fpe_env)==0) {      if (setjmp(g_fpe_env)==0) {
# Line 1118  void XASCWV( char *msg, /* pointer to st Line 1130  void XASCWV( char *msg, /* pointer to st
1130               double *r1,               double *r1,
1131               double *r2               double *r2
1132  ){  ){
1133        static double r1last;
1134    
1135        asc_assert(*level!=2); // LSODE doesn't give level 2 in our version.
1136    
1137      switch(*nerr){      switch(*nerr){
1138            case 52:
1139                if(*nr==2){
1140                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Illegal t = %f, not in range (t - hu,t) = (%f,%f)", r1last, *r1, *r2);
1141                    return;
1142                }else if(*nr==1){
1143                    r1last = *r1;
1144                    return;
1145                } break;
1146          case 204:          case 204:
1147              if(*nr==0)return;              if(*nr==2){
1148              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);
1149              break;                  return;
1150                } break;
1151          case 205:          case 205:
1152              if(*nr==0)return;              if(*nr==2){
1153              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Corrector convergence test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Corrector convergence test failed repeatedly or with abs(h)=hmin.\nt=%f and step size h=%f",*r1,*r2);
1154              break;                  return;
1155                } break;
1156            case 27:
1157                if(*nr==1 && *ni==1){
1158                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Trouble with INTDY: itask = %d, tout = %f", *i1, *r1);
1159                    return;
1160                } break;
1161        }
1162    
1163          default:      ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
             ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);  
1164    
1165              /* note that %.*s means that a string length (integer) and string pointer are being required */      /* note that %.*s means that a string length (integer) and string pointer are being required */
1166              FPRINTF(stderr,"LSODE error: %.*s",*nmes,msg);      FPRINTF(stderr,"LSODE error: (%d) %.*s",*nerr,*nmes,msg);
1167              if (*ni == 1) {      if (*ni == 1) {
1168              FPRINTF(stderr,"\nwhere i1 = %d\n",*i1);      FPRINTF(stderr,"\nwhere i1 = %d",*i1);
             }  
             if (*ni == 2) {  
             FPRINTF(stderr,"\nwhere i1 = %d, i2 = %d",*i1,*i2);  
             }  
             if (*nr == 1) {  
             FPRINTF(stderr,"\nwhere r1 = %.13g", *r1);  
             }  
             if (*nr == 2) {  
             FPRINTF(stderr,"\nwhere r1 = %.13g, r2 = %.13g", *r1,*r2);  
             }  
             error_reporter_end_flush();  
1169      }      }
1170        if (*ni == 2) {
1171      if (*level != 2) {      FPRINTF(stderr,"\nwhere i1 = %d, i2 = %d",*i1,*i2);
1172          return;      }
1173        if (*nr == 1) {
1174        FPRINTF(stderr,"\nwhere r1 = %.13g", *r1);
1175        }
1176        if (*nr == 2) {
1177        FPRINTF(stderr,"\nwhere r1 = %.13g, r2 = %.13g", *r1,*r2);
1178      }      }
   
     /* NOT reached. lsode does NOT make level 2 calls in our version. */  
1179      error_reporter_end_flush();      error_reporter_end_flush();
     Asc_Panic(3,"xascwv", "LSODE really really confused");  
1180  }  }
Legend:
Removed from v.965  
changed lines
  Added in v.966
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22