/[ascend]/trunk/base/generic/integrator/ida.c
ViewVC logotype

Diff of /trunk/base/generic/integrator/ida.c

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

revision 1262 by johnpye, Sat Jan 27 10:50:40 2007 UTC revision 1263 by johnpye, Sun Jan 28 04:03:01 2007 UTC
# Line 153  typedef struct IntegratorIdaPrecDJStruct Line 153  typedef struct IntegratorIdaPrecDJStruct
153      Hold all the function pointers associated with a particular preconditioner      Hold all the function pointers associated with a particular preconditioner
154      We don't need to store the 'pfree' function here as it is allocated to the enginedata struct      We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
155      by the pcreate function (ensures that corresponding 'free' and 'create' are always used)      by the pcreate function (ensures that corresponding 'free' and 'create' are always used)
156        
157      @note IDA uses a different convention for function pointer types, so no '*'.      @note IDA uses a different convention for function pointer types, so no '*'.
158  */  */
159  typedef struct IntegratorIdaPrecStruct{  typedef struct IntegratorIdaPrecStruct{
# Line 207  void integrator_ida_write_stats(Integrat Line 207  void integrator_ida_write_stats(Integrat
207  void integrator_ida_write_incidence(IntegratorSystem *blsys);  void integrator_ida_write_incidence(IntegratorSystem *blsys);
208    
209  /*------  /*------
210    Jacobi preconditioner -- experimental    Jacobi preconditioner -- experimental
211  */  */
212    
213  int integrator_ida_psetup_jacobi(realtype tt,  int integrator_ida_psetup_jacobi(realtype tt,
# Line 745  int integrator_ida_solve( Line 745  int integrator_ida_solve(
745          }          }
746          if (setjmp(g_fpe_env)==0) {          if (setjmp(g_fpe_env)==0) {
747  #endif  #endif
748            
749  # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3  # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
750          flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3  */          flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3  */
751  # else  # else
# Line 870  int integrator_ida_solve( Line 870  int integrator_ida_solve(
870    RESIDUALS AND JACOBIAN    RESIDUALS AND JACOBIAN
871  */  */
872    
873    #if 0
874  typedef void (SignalHandlerFn)(int);  typedef void (SignalHandlerFn)(int);
875  SignalHandlerFn integrator_ida_sig;  SignalHandlerFn integrator_ida_sig;
876  SignalHandlerFn *integrator_ida_sig_old;  SignalHandlerFn *integrator_ida_sig_old;
# Line 899  void integrator_ida_sig(int sig){ Line 900  void integrator_ida_sig(int sig){
900      CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);      CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);
901      longjmp(integrator_ida_jmp_buf,sig);      longjmp(integrator_ida_jmp_buf,sig);
902  }  }
903    #endif
904    
905  /**  /**
906      Function to evaluate system residuals, in the form required for IDA.      Function to evaluate system residuals, in the form required for IDA.
# Line 966  int integrator_ida_fex(realtype tt, N_Ve Line 968  int integrator_ida_fex(realtype tt, N_Ve
968      if (SETJMP(g_fpe_env)==0) {      if (SETJMP(g_fpe_env)==0) {
969  #endif  #endif
970    
971        
972  /*  /*
973      integrator_ida_sig_old = signal(SIGFPE,&integrator_ida_sig);      integrator_ida_sig_old = signal(SIGFPE,&integrator_ida_sig);
974      if(fegetenv(&integrator_ida_fenv_old))ASC_PANIC("fenv problem");      if(fegetenv(&integrator_ida_fenv_old))ASC_PANIC("fenv problem");
# Line 999  int integrator_ida_fex(realtype tt, N_Ve Line 1001  int integrator_ida_fex(realtype tt, N_Ve
1001          break;          break;
1002      default:      default:
1003          ERROR_REPORTER_HERE(ASC_PROG_ERR,"What exception?");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"What exception?");
1004          is_error=1;              is_error=1;
1005      }      }
1006      if(signal(SIGFPE,integrator_ida_sig_old)!=&integrator_ida_sig)ASC_PANIC("signal problem");      if(signal(SIGFPE,integrator_ida_sig_old)!=&integrator_ida_sig)ASC_PANIC("signal problem");
1007      if(fesetenv(&integrator_ida_fenv_old))ASC_PANIC("fenv problem");      if(fesetenv(&integrator_ida_fenv_old))ASC_PANIC("fenv problem");
# Line 1348  int integrator_ida_jvex(realtype tt, N_V Line 1350  int integrator_ida_jvex(realtype tt, N_V
1350                  if(variables[j] == blsys->x) continue;                  if(variables[j] == blsys->x) continue;
1351  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
1352                  CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));                  CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1353  #endif                #endif
1354                  if(var_deriv(variables[j])){                  if(var_deriv(variables[j])){
1355  #define DIFFINDEX integrator_ida_diffindex(blsys,variables[j])  #define DIFFINDEX integrator_ida_diffindex(blsys,variables[j])
1356  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
# Line 1592  enum integrator_ida_write_jac_enum{ Line 1594  enum integrator_ida_write_jac_enum{
1594  */  */
1595  void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){  void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){
1596      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
1597      MM_typecode matcode;                              MM_typecode matcode;
1598      int nnz, rhomax;      int nnz, rhomax;
1599      double *derivatives;      double *derivatives;
1600      struct var_variable **variables;      struct var_variable **variables;
# Line 1623  void integrator_ida_write_jacobian(Integ Line 1625  void integrator_ida_write_jacobian(Integ
1625      mm_set_matrix(&matcode);      mm_set_matrix(&matcode);
1626      mm_set_coordinate(&matcode);      mm_set_coordinate(&matcode);
1627      mm_set_real(&matcode);      mm_set_real(&matcode);
1628      mm_write_banner(f, matcode);      mm_write_banner(f, matcode);
1629      mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);          mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1630    
1631      variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);      variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1632      derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);      derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
# Line 1656  void integrator_ida_write_jacobian(Integ Line 1658  void integrator_ida_write_jacobian(Integ
1658                  if(!var_deriv(variables[j])){                  if(!var_deriv(variables[j])){
1659                      fprintf(f, "%d %d %10.3g\n", i, var_sindex(variables[j]), derivatives[j]);                      fprintf(f, "%d %d %10.3g\n", i, var_sindex(variables[j]), derivatives[j]);
1660                  }                  }
1661              }                                            }
1662          }          }
1663      }      }
1664      ASC_FREE(variables);      ASC_FREE(variables);
# Line 1700  void integrator_ida_write_incidence(Inte Line 1702  void integrator_ida_write_incidence(Inte
1702              break;              break;
1703          }          }
1704    
1705          fprintf(stderr,"%3d:%-15s:",i,relname);              fprintf(stderr,"%3d:%-15s:",i,relname);
1706          ASC_FREE(relname);          ASC_FREE(relname);
1707    
1708          for(j=0; j<count; ++j){          for(j=0; j<count; ++j){
# Line 1792  int integrator_ida_debug(const Integrato Line 1794  int integrator_ida_debug(const Integrato
1794      /* let's write out the relations too */      /* let's write out the relations too */
1795      rlist = slv_get_solvers_rel_list(sys->system);      rlist = slv_get_solvers_rel_list(sys->system);
1796      rlen = slv_get_num_solvers_rels(sys->system);      rlen = slv_get_num_solvers_rels(sys->system);
1797        
1798      fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);      fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
1799      fprintf(fp,"index\tname\n");      fprintf(fp,"index\tname\n");
1800      fprintf(fp,"-----\t----\n");      fprintf(fp,"-----\t----\n");
# Line 1810  int integrator_ida_debug(const Integrato Line 1812  int integrator_ida_debug(const Integrato
1812          return 340;          return 340;
1813      }      }
1814      fprintf(fp,"\n");      fprintf(fp,"\n");
1815        
1816      /* and lets write block debug output */      /* and lets write block debug output */
1817      system_block_debug(sys->system, fp);      system_block_debug(sys->system, fp);
1818            
1819      return 0; /* success */      return 0; /* success */
1820  }  }
1821    
Legend:
Removed from v.1262  
changed lines
  Added in v.1263
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22