/[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 1048 by johnpye, Sun Dec 24 01:33:59 2006 UTC revision 1049 by johnpye, Fri Jan 5 13:45:13 2007 UTC
# Line 680  static void LSODE_FEX( int *n_eq ,double Line 680  static void LSODE_FEX( int *n_eq ,double
680    
681    /*  slv_parameters_t parameters; pity lsode doesn't allow error returns */    /*  slv_parameters_t parameters; pity lsode doesn't allow error returns */
682    /* int i; */    /* int i; */
683    unsigned long ok;    unsigned long res;
684    
685  #if DOTIME  #if DOTIME
686    double time1,time2;    double time1,time2;
# Line 726  static void LSODE_FEX( int *n_eq ,double Line 726  static void LSODE_FEX( int *n_eq ,double
726    slv_solve(l_lsode_blsys->system);    slv_solve(l_lsode_blsys->system);
727    slv_get_status(l_lsode_blsys->system, &status);    slv_get_status(l_lsode_blsys->system, &status);
728    /* pass the solver status to the integrator */    /* pass the solver status to the integrator */
729    ok = integrator_checkstatus(status);    res = integrator_checkstatus(status);
730    
731  #if DOTIME  #if DOTIME
732    time2 = tm_cpu_time() - time2;    time2 = tm_cpu_time() - time2;
733  #endif  #endif
734    
735    if (!ok) {    if(res){
736      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for derivatives");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for derivatives (%d)",res);
737      /*      /*
738      ERROR_REPORTER_START_HERE(ASC_PROG_ERR);      ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
739      FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");      FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");
# Line 743  static void LSODE_FEX( int *n_eq ,double Line 743  static void LSODE_FEX( int *n_eq ,double
743      error_reporter_end_flush();      error_reporter_end_flush();
744      */      */
745      lsodesys.status = lsode_nok;      lsodesys.status = lsode_nok;
746    } else {    }else{
747      lsodesys.status = lsode_ok;      lsodesys.status = lsode_ok;
748    }    }
749    integrator_get_ydot(l_lsode_blsys, ydot);    integrator_get_ydot(l_lsode_blsys, ydot);
# Line 820  static void LSODE_JEX(int *neq ,double * Line 820  static void LSODE_JEX(int *neq ,double *
820  /**  /**
821      The public function: here we do the actual integration, I guess.      The public function: here we do the actual integration, I guess.
822    
823      Return 1 on success      Return 0 on success
824  */  */
825  int integrator_lsode_solve(IntegratorSystem *blsys  int integrator_lsode_solve(IntegratorSystem *blsys
826          , unsigned long start_index, unsigned long finish_index          , unsigned long start_index, unsigned long finish_index
# Line 866  int integrator_lsode_solve(IntegratorSys Line 866  int integrator_lsode_solve(IntegratorSys
866       We handle any linsol/linsolqr based solver. */       We handle any linsol/linsolqr based solver. */
867    if (strcmp(slv_solver_name(slv_get_selected_solver(blsys->system)),"QRSlv") != 0) {    if (strcmp(slv_solver_name(slv_get_selected_solver(blsys->system)),"QRSlv") != 0) {
868      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv must be selected before integration.");      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv must be selected before integration.");
869      return 0;      return 1;
870    }    }
871    
872    slv_get_status(l_lsode_blsys->system, &status);    slv_get_status(l_lsode_blsys->system, &status);
# Line 874  int integrator_lsode_solve(IntegratorSys Line 874  int integrator_lsode_solve(IntegratorSys
874    if (status.struct_singular) {    if (status.struct_singular) {
875      ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system is structurally singular.");      ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system is structurally singular.");
876      lsodesys.status = lsode_nok;      lsodesys.status = lsode_nok;
877      return 0;      return 2;
878    }    }
879    
880  #if defined(STATIC_LSOD) || defined (DYNAMIC_LSOD)  #if defined(STATIC_LSOD) || defined (DYNAMIC_LSOD)
# Line 888  int integrator_lsode_solve(IntegratorSys Line 888  int integrator_lsode_solve(IntegratorSys
888    if (nsamples <2) {    if (nsamples <2) {
889      ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined.");      ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined.");
890      lsodesys.status = lsode_nok;      lsodesys.status = lsode_nok;
891      return 0;      return 3;
892    }    }
893    neq = blsys->n_y;    neq = blsys->n_y;
894    nobs = blsys->n_obs;    nobs = blsys->n_obs;
# Line 911  int integrator_lsode_solve(IntegratorSys Line 911  int integrator_lsode_solve(IntegratorSys
911      lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);      lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
912      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for lsode.");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for lsode.");
913      lsodesys.status = lsode_nok;      lsodesys.status = lsode_nok;
914      return 0;      return 4;
915    }    }
916    
917    /*    /*
# Line 929  int integrator_lsode_solve(IntegratorSys Line 929  int integrator_lsode_solve(IntegratorSys
929    
930    if(x[0] > integrator_getsample(blsys, 2)){    if(x[0] > integrator_getsample(blsys, 2)){
931      ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value");      ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value");
932      return 0;      return 5;
933    }    }
934    
935    /* put the values from derivative system into the record */    /* put the values from derivative system into the record */
# Line 998  int integrator_lsode_solve(IntegratorSys Line 998  int integrator_lsode_solve(IntegratorSys
998        if (obs_out!=NULL) {        if (obs_out!=NULL) {
999          fclose(obs_out);          fclose(obs_out);
1000        }        }
1001        return 0;        return 6;
1002      }      }
1003  # endif /* NO_SIGNAL_TRAPS */  # endif /* NO_SIGNAL_TRAPS */
1004    
# Line 1019  int integrator_lsode_solve(IntegratorSys Line 1019  int integrator_lsode_solve(IntegratorSys
1019    
1020        lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);        lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
1021        integrator_output_close(blsys);        integrator_output_close(blsys);
1022        return 0;        return 7;
1023      }      }
1024    
1025      if (lsodesys.status==lsode_nok) {      if (lsodesys.status==lsode_nok) {
# Line 1028  int integrator_lsode_solve(IntegratorSys Line 1028  int integrator_lsode_solve(IntegratorSys
1028        lsodesys.status = lsode_ok;       /* clean up before we go */        lsodesys.status = lsode_ok;       /* clean up before we go */
1029        lsodesys.lastcall = lsode_none;        lsodesys.lastcall = lsode_none;
1030        integrator_output_close(blsys);        integrator_output_close(blsys);
1031        return 0;        return 8;
1032      }      }
1033    
1034      integrator_setsample(blsys, index+1, x[0]);      integrator_setsample(blsys, index+1, x[0]);
# Line 1045  int integrator_lsode_solve(IntegratorSys Line 1045  int integrator_lsode_solve(IntegratorSys
1045          lsodesys.status = lsode_ok;          lsodesys.status = lsode_ok;
1046          lsodesys.lastcall = lsode_none;          lsodesys.lastcall = lsode_none;
1047          integrator_output_close(blsys);          integrator_output_close(blsys);
1048          return 0;          return 9;
1049      }      }
1050    
1051      if (nobs > 0) {      if (nobs > 0) {
# Line 1073  int integrator_lsode_solve(IntegratorSys Line 1073  int integrator_lsode_solve(IntegratorSys
1073          lsodesys.status = lsode_ok;               /* clean up before we go */          lsodesys.status = lsode_ok;               /* clean up before we go */
1074          lsodesys.lastcall = lsode_none;          lsodesys.lastcall = lsode_none;
1075          integrator_output_close(blsys);          integrator_output_close(blsys);
1076          return 0;          return 10;
1077        }        }
1078  # endif /* NO_SIGNAL_TRAPS */  # endif /* NO_SIGNAL_TRAPS */
1079      }      }
# Line 1099  int integrator_lsode_solve(IntegratorSys Line 1099  int integrator_lsode_solve(IntegratorSys
1099    integrator_output_close(blsys);    integrator_output_close(blsys);
1100    
1101    CONSOLE_DEBUG("--- LSODE done ---");    CONSOLE_DEBUG("--- LSODE done ---");
1102    return 1;    return 0; /* success */
1103    
1104  #else /* STATIC_LSOD || DYNAMIC_LSOD */  #else /* STATIC_LSOD || DYNAMIC_LSOD */
1105    
1106    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration will not be performed. LSODE binary not available.");    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration will not be performed. LSODE binary not available.");
1107    lsodesys.status = lsode_nok;    lsodesys.status = lsode_nok;
1108    return 0;    return 11;
1109    
1110  #endif  #endif
1111  }  }
# Line 1153  void XASCWV( char *msg, /* pointer to st Line 1153  void XASCWV( char *msg, /* pointer to st
1153                  return;                  return;
1154              } break;              } break;
1155          case 204:          case 204:
1156                if(*nr==0 && *ni==0)return;
1157              if(*nr==2){              if(*nr==2){
1158                  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);
1159                  return;                  return;
Legend:
Removed from v.1048  
changed lines
  Added in v.1049
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22