/[ascend]/trunk/solvers/ida/ida.c
ViewVC logotype

Diff of /trunk/solvers/ida/ida.c

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

revision 2368 by jpye, Mon Jan 31 06:41:54 2011 UTC revision 2369 by jpye, Mon Jan 31 09:00:44 2011 UTC
# Line 89  Line 89 
89  #include "idalinear.h"  #include "idalinear.h"
90  #include "idaanalyse.h"  #include "idaanalyse.h"
91  #include "ida_impl.h"  #include "ida_impl.h"
92    #include "idaprec.h"
93    #include "idacalc.h"
94    
95  /*  /*
96      for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2      for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
# Line 158  extern ASC_EXPORT int ida_register(void) Line 160  extern ASC_EXPORT int ida_register(void)
160    FORWARD DECLS    FORWARD DECLS
161  */  */
162    
 /* forward dec needed for IntegratorIdaPrecFreeFn */  
 struct IntegratorIdaDataStruct;  
   
 /* functions for allocating storage for and freeing preconditioner data */  
 typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *integ);  
 typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);  
   
   
 /**  
     Struct containing any stuff that IDA needs that doesn't fit into the  
     common IntegratorSystem struct.  
 */  
 typedef struct IntegratorIdaDataStruct{  
   
     struct rel_relation **rellist;   /**< NULL terminated list of ACTIVE rels */  
     int nrels; /* number of ACTIVE rels */  
   
     struct bnd_boundary **bndlist;   /**< NULL-terminated list of boundaries, for use in the root-finding  code */  
     int nbnds; /* number of boundaries */  
   
     int safeeval;                    /**< whether to pass the 'safe' flag to relman_eval */  
     var_filter_t vfilter;  
     rel_filter_t rfilter;            /**< Used to filter relations from solver's rellist (@TODO needs work) */  
     void *precdata;                  /**< For use by the preconditioner */  
     IntegratorIdaPrecFreeFn *pfree;  /**< Store instructions here on how to free precdata */  
   
 } IntegratorIdaData;  
   
   
 typedef struct IntegratorIdaPrecDJStruct{  
     N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */  
 } IntegratorIdaPrecDataJacobi;  
   
 typedef struct IntegratorIdaPrecDJFStruct{  
     linsolqr_system_t L;  
 } IntegratorIdaPrecDataJacobian;  
   
 /**  
     Hold all the function pointers associated with a particular preconditioner  
     We don't need to store the 'pfree' function here as it is allocated to the enginedata struct  
     by the pcreate function (ensures that corresponding 'free' and 'create' are always used)  
   
     @note IDA uses a different convention for function pointer types, so no '*'.  
 */  
 typedef struct IntegratorIdaPrecStruct{  
     IntegratorIdaPrecCreateFn *pcreate;  
     IDASpilsPrecSetupFn psetup;  
     IDASpilsPrecSolveFn psolve;  
 } IntegratorIdaPrec;  
   
 /* residual function forward declaration */  
 static int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);  
   
 static int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr  
         , N_Vector v, N_Vector Jv, realtype c_j  
         , void *jac_data, N_Vector tmp1, N_Vector tmp2  
 );  
   
163  /* error handler forward declaration */  /* error handler forward declaration */
164  static void integrator_ida_error(int error_code  static void integrator_ida_error(int error_code
165          , const char *module, const char *function          , const char *module, const char *function
166          , char *msg, void *eh_data          , char *msg, void *eh_data
167  );  );
168    
 /* dense jacobian evaluation for IDADense dense direct linear solver */  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
 static int integrator_ida_djex(int Neq, realtype tt, realtype c_j  
         , N_Vector yy, N_Vector yp, N_Vector rr  
         , IDA_MTX_T Jac, void *jac_data  
         , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3  
 );  
 #else  
 static int integrator_ida_djex(long int Neq, realtype tt  
         , N_Vector yy, N_Vector yp, N_Vector rr  
         , realtype c_j, void *jac_data, IDA_MTX_T Jac  
         , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3  
 );  
 #endif  
   
 /* sparse jacobian evaluation for ASCEND's sparse direct solver */  
 static IntegratorSparseJacFn integrator_ida_sjex;  
   
 /* boundary-detection function */  
 static int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);  
   
169  typedef struct IntegratorIdaStatsStruct{  typedef struct IntegratorIdaStatsStruct{
170      long nsteps;      long nsteps;
171      long nrevals;      long nrevals;
# Line 261  static void integrator_dae_show_var(Inte Line 184  static void integrator_dae_show_var(Inte
184    
185  static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);  static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
186  static void integrator_ida_write_stats(IntegratorIdaStats *stats);  static void integrator_ida_write_stats(IntegratorIdaStats *stats);
 static void integrator_ida_write_incidence(IntegratorSystem *integ);  
   
 /*------  
   Full jacobian preconditioner -- experimental  
 */  
   
 static int integrator_ida_psetup_jacobian(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          realtype c_j, void *prec_data,  
          N_Vector tmp1, N_Vector tmp2,  
          N_Vector tmp3  
 );  
   
 static int integrator_ida_psolve_jacobian(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          N_Vector rvec, N_Vector zvec,  
          realtype c_j, realtype delta, void *prec_data,  
          N_Vector tmp  
 );  
   
 static void integrator_ida_pcreate_jacobian(IntegratorSystem *integ);  
   
 static void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata);  
   
 static const IntegratorIdaPrec prec_jacobian = {  
     integrator_ida_pcreate_jacobian  
     , integrator_ida_psetup_jacobian  
     , integrator_ida_psolve_jacobian  
 };  
   
 /*------  
   Jacobi preconditioner -- experimental  
 */  
   
 static int integrator_ida_psetup_jacobi(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          realtype c_j, void *prec_data,  
          N_Vector tmp1, N_Vector tmp2,  
          N_Vector tmp3  
 );  
   
 static int integrator_ida_psolve_jacobi(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          N_Vector rvec, N_Vector zvec,  
          realtype c_j, realtype delta, void *prec_data,  
          N_Vector tmp  
 );  
187    
 static void integrator_ida_pcreate_jacobi(IntegratorSystem *integ);  
   
 static void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);  
   
 static const IntegratorIdaPrec prec_jacobi = {  
     integrator_ida_pcreate_jacobi  
     , integrator_ida_psetup_jacobi  
     , integrator_ida_psolve_jacobi  
 };  
188    
189  /*-------------------------------------------------------------  /*-------------------------------------------------------------
190    SETUP/TEARDOWN ROUTINES    SETUP/TEARDOWN ROUTINES
# Line 364  static void integrator_ida_free(void *en Line 231  static void integrator_ida_free(void *en
231  #endif  #endif
232  }  }
233    
234  static IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *integ){  IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *integ){
235      IntegratorIdaData *d;      IntegratorIdaData *d;
236      assert(integ!=NULL);      assert(integ!=NULL);
237      assert(integ->enginedata!=NULL);      assert(integ->enginedata!=NULL);
# Line 1156  static int integrator_ida_solve( Line 1023  static int integrator_ida_solve(
1023      return 0;      return 0;
1024  }  }
1025    
 /*--------------------------------------------------  
   RESIDUALS AND JACOBIAN AND IDAROOTFN  
 */  
   
 #if 0  
 typedef void (SignalHandlerFn)(int);  
 SignalHandlerFn integrator_ida_sig;  
 SignalHandlerFn *integrator_ida_sig_old;  
 jmp_buf integrator_ida_jmp_buf;  
 fenv_t integrator_ida_fenv_old;  
   
   
 void integrator_ida_write_feinfo(){  
     int f;  
     f = fegetexcept();  
     CONSOLE_DEBUG("Locating nature of exception...");  
     if(f & FE_DIVBYZERO)ERROR_REPORTER_HERE(ASC_PROG_ERR,"DIV BY ZERO");  
     if(f & FE_INEXACT)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INEXACT");  
     if(f & FE_INVALID)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INVALID");  
     if(f & FE_OVERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"OVERFLOW");  
     if(f & FE_UNDERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"UNDERFLOW");  
     if(f==0)ERROR_REPORTER_HERE(ASC_PROG_ERR,"FLAGS ARE CLEAR?!?");  
 }  
   
 void integrator_ida_sig(int sig){  
     /* the wrong signal: rethrow to the default handler */  
     if(sig!=SIGFPE){  
         signal(SIGFPE,SIG_DFL);  
         raise(sig);  
     }  
     integrator_ida_write_feinfo();  
     CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);  
     longjmp(integrator_ida_jmp_buf,sig);  
 }  
 #endif  
   
 /**  
     Function to evaluate system residuals, in the form required for IDA.  
   
     Given tt, yy and yp, we need to evaluate and return rr.  
   
     @param tt current value of indep variable (time)  
     @param yy current values of dependent variable vector  
     @param yp current values of derivatives of dependent variables  
     @param rr the output residual vector (is we're returning data to)  
     @param res_data pointer to our stuff (integ in this case).  
   
     @return 0 on success, positive on recoverable error, and  
         negative on unrecoverable error.  
 */  
 static int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){  
     IntegratorSystem *integ;  
     IntegratorIdaData *enginedata;  
     int i, calc_ok, is_error;  
     struct rel_relation** relptr;  
     double resid;  
     char *relname;  
 #ifdef FEX_DEBUG  
     char *varname;  
     char diffname[30];  
 #endif  
   
     integ = (IntegratorSystem *)res_data;  
     enginedata = integrator_ida_enginedata(integ);  
   
 #ifdef FEX_DEBUG  
     /* fprintf(stderr,"\n\n"); */  
     CONSOLE_DEBUG("EVALUTE RESIDUALS...");  
 #endif  
   
     if(NV_LENGTH_S(rr)!=enginedata->nrels){  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");  
         return -1; /* unrecoverable */  
     }  
   
     /* pass the values of everything back to the compiler */  
     integrator_set_t(integ, (double)tt);  
     integrator_set_y(integ, NV_DATA_S(yy));  
     integrator_set_ydot(integ, NV_DATA_S(yp));  
   
     /* perform bounds checking on all variables */  
     if(slv_check_bounds(integ->system, 0, -1, NULL)){  
         /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */  
         return 1;  
     }  
   
     /* evaluate each residual in the rellist */  
     is_error = 0;  
     relptr = enginedata->rellist;  
   
   
 #ifdef ASC_SIGNAL_TRAPS  
     if(enginedata->safeeval){  
         Asc_SignalHandlerPush(SIGFPE,SIG_IGN);  
     }else{  
 # ifdef FEX_DEBUG  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");  
 # endif  
         Asc_SignalHandlerPushDefault(SIGFPE);  
     }  
   
     if (SETJMP(g_fpe_env)==0) {  
 #endif  
   
   
     for(i=0, relptr = enginedata->rellist;  
                 i< enginedata->nrels && relptr != NULL;  
                 ++i, ++relptr  
     ){  
         resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);  
   
         NV_Ith_S(rr,i) = resid;  
         if(!calc_ok){  
             relname = rel_make_name(integ->system, *relptr);  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);  
             ASC_FREE(relname);  
             /* presumable some output already made? */  
             is_error = 1;  
         }/*else{  
             CONSOLE_DEBUG("Calc OK");  
         }*/  
     }  
   
     if(!is_error){  
         for(i=0;i< enginedata->nrels; ++i){  
             if(isnan(NV_Ith_S(rr,i))){  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);  
                 is_error=1;  
             }  
         }  
 #ifdef FEX_DEBUG  
         if(!is_error){  
             CONSOLE_DEBUG("No NAN detected");  
         }  
 #endif  
     }  
   
 #ifdef ASC_SIGNAL_TRAPS  
     }else{  
         relname = rel_make_name(integ->system, *relptr);  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);  
         ASC_FREE(relname);  
         is_error = 1;  
     }  
   
     if(enginedata->safeeval){  
         Asc_SignalHandlerPop(SIGFPE,SIG_IGN);  
     }else{  
         Asc_SignalHandlerPopDefault(SIGFPE);  
     }  
 #endif  
   
   
 #ifdef FEX_DEBUG  
     /* output residuals to console */  
     CONSOLE_DEBUG("RESIDUAL OUTPUT");  
     fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");  
     for(i=0; i<integ->n_y; ++i){  
         varname = var_make_name(integ->system,integ->y[i]);  
         fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));  
         if(integ->ydot[i]){  
             varname = var_make_name(integ->system,integ->ydot[i]);  
             fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));  
         }else{  
             snprintf(diffname,99,"diff(%s)",varname);  
             fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));  
         }  
         ASC_FREE(varname);  
         relname = rel_make_name(integ->system,enginedata->rellist[i]);  
         fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);  
     }  
 #endif  
   
     if(is_error){  
         return 1;  
     }  
   
 #ifdef FEX_DEBUG  
     CONSOLE_DEBUG("RESIDUAL OK");  
 #endif  
     return 0;  
 }  
   
 /**  
     Dense Jacobian evaluation. Only suitable for small problems!  
     Has been seen working for problems up to around 2000 vars, FWIW.  
 */  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
 static int integrator_ida_djex(int Neq, realtype tt, realtype c_j  
         , N_Vector yy, N_Vector yp, N_Vector rr  
         , IDA_MTX_T Jac, void *jac_data  
         , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3  
 ){  
 #else  
 static int integrator_ida_djex(long int Neq, realtype tt  
         , N_Vector yy, N_Vector yp, N_Vector rr  
         , realtype c_j, void *jac_data, IDA_MTX_T Jac  
         , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3  
 ){  
 #endif  
     IntegratorSystem *integ;  
     IntegratorIdaData *enginedata;  
     char *relname;  
 #ifdef DJEX_DEBUG  
     struct var_variable **varlist;  
     char *varname;  
 #endif  
     struct rel_relation **relptr;  
     int i;  
     double *derivatives;  
     struct var_variable **variables;  
     int count, j;  
     int status, is_error = 0;  
   
     integ = (IntegratorSystem *)jac_data;  
     enginedata = integrator_ida_enginedata(integ);  
   
     /* allocate space for returns from relman_diff3 */  
     /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */  
     variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);  
     derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);  
   
     /* pass the values of everything back to the compiler */  
     integrator_set_t(integ, (double)tt);  
     integrator_set_y(integ, NV_DATA_S(yy));  
     integrator_set_ydot(integ, NV_DATA_S(yp));  
   
     /* perform bounds checking on all variables */  
     if(slv_check_bounds(integ->system, 0, -1, NULL)){  
         /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */  
         return 1;  
     }  
   
 #ifdef DJEX_DEBUG  
     varlist = slv_get_solvers_var_list(integ->system);  
   
     /* print vars */  
     for(i=0; i < integ->n_y; ++i){  
         varname = var_make_name(integ->system, integ->y[i]);  
         CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));  
         asc_assert(NV_Ith_S(yy,i) == var_value(integ->y[i]));  
         ASC_FREE(varname);  
     }  
   
     /* print derivatives */  
     for(i=0; i < integ->n_y; ++i){  
         if(integ->ydot[i]){  
             varname = var_make_name(integ->system, integ->ydot[i]);  
             CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(integ->ydot[i]));  
             ASC_FREE(varname);  
         }else{  
             varname = var_make_name(integ->system, integ->y[i]);  
             CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));  
             ASC_FREE(varname);  
         }  
     }  
   
     /* print step size */  
     CONSOLE_DEBUG("<c_j> = %g",c_j);  
 #endif  
   
     /* build up the dense jacobian matrix... */  
     status = 0;  
     for(i=0, relptr = enginedata->rellist;  
             i< enginedata->nrels && relptr != NULL;  
             ++i, ++relptr  
     ){  
         /* get derivatives for this particular relation */  
         status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);  
   
         if(status){  
             relname = rel_make_name(integ->system, *relptr);  
             CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);  
             ASC_FREE(relname);  
             is_error = 1;  
             break;  
         }  
   
         /* output what's going on here ... */  
 #ifdef DJEX_DEBUG  
         relname = rel_make_name(integ->system, *relptr);  
         fprintf(stderr,"%d: '%s': ",i,relname);  
         for(j=0;j<count;++j){  
             varname = var_make_name(integ->system, variables[j]);  
             if(var_deriv(variables[j])){  
                 fprintf(stderr,"  '%s'=",varname);  
                 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(integ,variables[j]));  
             }else{  
                 fprintf(stderr,"  '%s'=y[%d]",varname,var_sindex(variables[j]));  
             }  
             ASC_FREE(varname);  
         }  
         /* relname is freed further down */  
         fprintf(stderr,"\n");  
 #endif  
   
         /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */  
         for(j=0; j < count; ++j){  
 #ifdef DJEX_DEBUG  
             varname = var_make_name(integ->system,variables[j]);  
             fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);  
             ASC_FREE(varname);  
 #endif  
             if(!var_deriv(variables[j])){  
 #ifdef DJEX_DEBUG  
                 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);  
                 asc_assert(var_sindex(variables[j]) >= 0);  
                 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);  
 #endif  
                 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];  
             }else{  
                 DENSE_ELEM(Jac,i,integrator_ida_diffindex(integ,variables[j])) += derivatives[j] * c_j;  
 #ifdef DJEX_DEBUG  
                 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);  
 #endif  
             }  
         }  
     }  
   
 #ifdef DJEX_DEBUG  
     ASC_FREE(relname);  
     CONSOLE_DEBUG("PRINTING JAC");  
     fprintf(stderr,"\t");  
     for(j=0; j < integ->n_y; ++j){  
         if(j)fprintf(stderr,"\t");  
         varname = var_make_name(integ->system,integ->y[j]);  
         fprintf(stderr,"%11s",varname);  
         ASC_FREE(varname);  
     }  
     fprintf(stderr,"\n");  
     for(i=0; i < enginedata->nrels; ++i){  
         relname = rel_make_name(integ->system, enginedata->rellist[i]);  
         fprintf(stderr,"%s\t",relname);  
         ASC_FREE(relname);  
   
         for(j=0; j < integ->n_y; ++j){  
             if(j)fprintf(stderr,"\t");  
             fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));  
         }  
         fprintf(stderr,"\n");  
     }  
 #endif  
   
     /* test for NANs */  
     if(!is_error){  
         for(i=0;i< enginedata->nrels; ++i){  
             for(j=0;j<integ->n_y;++j){  
                 if(isnan(DENSE_ELEM(Jac,i,j))){  
                     ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);  
                     is_error=1;  
                 }  
             }  
         }  
 #ifdef DJEX_DEBUG  
         if(!is_error){  
             CONSOLE_DEBUG("No NAN detected");  
         }  
 #endif  
     }  
   
 /*  if(integrator_ida_check_diffindex(integ)){  
         is_error = 1;  
     }*/  
   
     if(is_error){  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");  
         return 1;  
     }  
   
 #ifdef DJEX_DEBUG  
     CONSOLE_DEBUG("DJEX RETURNING 0");  
     /* ASC_PANIC("Quitting"); */  
 #endif  
     return 0;  
 }  
   
 /**  
     Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)  
   
     Given tt, yy, yp, rr and v, we need to evaluate and return Jv.  
   
     @param tt current value of the independent variable (time, t)  
     @param yy current value of the dependent variable vector, y(t).  
     @param yp current value of y'(t).  
     @param rr current value of the residual vector F(t, y, y').  
     @param v  the vector by which the Jacobian must be multiplied to the right.  
     @param Jv the output vector computed  
     @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).  
     @param jac_data pointer to our stuff (integ in this case, passed into IDA via IDASp*SetJacTimesVecFn.)  
     @param tmp1 @see tmp2  
     @param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.  
     @return 0 on success  
 */  
 static int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr  
         , N_Vector v, N_Vector Jv, realtype c_j  
         , void *jac_data, N_Vector tmp1, N_Vector tmp2  
 ){  
     IntegratorSystem *integ;  
     IntegratorIdaData *enginedata;  
     int i, j, is_error=0;  
     struct rel_relation** relptr = 0;  
     char *relname;  
     int status;  
     double Jv_i;  
   
     struct var_variable **variables;  
     double *derivatives;  
     int count;  
     struct var_variable **varlist;  
 #ifdef JEX_DEBUG  
   
     CONSOLE_DEBUG("EVALUATING JACOBIAN...");  
 #endif  
   
     integ = (IntegratorSystem *)jac_data;  
     enginedata = integrator_ida_enginedata(integ);  
     varlist = slv_get_solvers_var_list(integ->system);  
   
     /* pass the values of everything back to the compiler */  
     integrator_set_t(integ, (double)tt);  
     integrator_set_y(integ, NV_DATA_S(yy));  
     integrator_set_ydot(integ, NV_DATA_S(yp));  
     /* no real use for residuals (rr) here, I don't think? */  
   
     /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */  
   
     i = NV_LENGTH_S(yy) * 2;  
 #ifdef JEX_DEBUG  
     CONSOLE_DEBUG("Allocating 'variables' with length %d",i);  
 #endif  
     variables = ASC_NEW_ARRAY(struct var_variable*, i);  
     derivatives = ASC_NEW_ARRAY(double, i);  
   
     /* evaluate the derivatives... */  
     /* J = dG_dy = dF_dy + alpha * dF_dyp */  
   
 #ifdef ASC_SIGNAL_TRAPS  
     Asc_SignalHandlerPushDefault(SIGFPE);  
     if (SETJMP(g_fpe_env)==0) {  
 #endif  
         for(i=0, relptr = enginedata->rellist;  
                 i< enginedata->nrels && relptr != NULL;  
                 ++i, ++relptr  
         ){  
             /* get derivatives for this particular relation */  
             status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);  
 #ifdef JEX_DEBUG  
             CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);  
 #endif  
   
             if(status){  
                 relname = rel_make_name(integ->system, *relptr);  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);  
                 ASC_FREE(relname);  
                 is_error = 1;  
                 break;  
             }  
   
             /*  
                 Now we have the derivatives wrt each alg/diff variable in the  
                 present equation. variables[] points into the varlist. need  
                 a mapping from the varlist to the y and ydot lists.  
             */  
   
             Jv_i = 0;  
             for(j=0; j < count; ++j){  
                 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], integ->n_y);  
                 varname = var_make_name(integ->system, enginedata->varlist[variables[j]]);  
                 if(varname){  
                     CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);  
                     ASC_FREE(varname);  
                 }else{  
                     CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);  
                 }  
                 */  
   
                 /* we don't calculate derivatives wrt indep var */  
                 asc_assert(variables[j]>=0);  
                 if(variables[j] == integ->x) continue;  
 #ifdef JEX_DEBUG  
                 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));  
 #endif  
                 if(var_deriv(variables[j])){  
 #define DIFFINDEX integrator_ida_diffindex(integ,variables[j])  
 #ifdef JEX_DEBUG  
                     fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i  
                         , derivatives[j] * NV_Ith_S(v,DIFFINDEX)  
                         , i, DIFFINDEX, derivatives[j]  
                         , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)  
                     );  
 #endif  
                     asc_assert(integ->ydot[DIFFINDEX]==variables[j]);  
                     Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;  
 #undef DIFFINDEX  
                 }else{  
 #define VARINDEX var_sindex(variables[j])  
 #ifdef JEX_DEBUG  
                     asc_assert(integ->y[VARINDEX]==variables[j]);  
                     fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"  
                         , i  
                         , derivatives[j] * NV_Ith_S(v,VARINDEX)  
                         , i, VARINDEX, derivatives[j]  
                         , VARINDEX, NV_Ith_S(v,VARINDEX)  
                     );  
 #endif  
                     Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);  
 #undef VARINDEX  
                 }  
             }  
   
             NV_Ith_S(Jv,i) = Jv_i;  
 #ifdef JEX_DEBUG  
             CONSOLE_DEBUG("rel = %p",*relptr);  
             relname = rel_make_name(integ->system, *relptr);  
             CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));  
             ASC_FREE(relname);  
             return 1;  
 #endif  
         }  
 #ifdef ASC_SIGNAL_TRAPS  
     }else{  
         relname = rel_make_name(integ->system, *relptr);  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);  
         ASC_FREE(relname);  
         is_error = 1;  
     }  
     Asc_SignalHandlerPopDefault(SIGFPE);  
 #endif  
   
     if(is_error){  
         CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");  
         return 1;  
     }  
     return 0;  
 }  
   
 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */  
 static int integrator_ida_sjex(long int Neq, realtype tt  
         , N_Vector yy, N_Vector yp, N_Vector rr  
         , realtype c_j, void *jac_data, mtx_matrix_t Jac  
         , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3  
 ){  
     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");  
     return -1;  
 }  
   
 /* root finding function */  
   
 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data){  
     IntegratorSystem *integ;  
     IntegratorIdaData *enginedata;  
     int i;  
 #ifdef ROOT_DEBUG  
     char *relname;  
 #endif  
   
     asc_assert(g_data!=NULL);  
     integ = (IntegratorSystem *)g_data;  
     enginedata = integrator_ida_enginedata(integ);  
   
     /* pass the values of everything back to the compiler */  
     integrator_set_t(integ, (double)tt);  
     integrator_set_y(integ, NV_DATA_S(yy));  
     integrator_set_ydot(integ, NV_DATA_S(yp));  
   
     asc_assert(gout!=NULL);  
   
 #ifdef ROOT_DEBUG  
     CONSOLE_DEBUG("t = %f",tt);  
 #endif  
   
     /* evaluate the residuals for each of the boundaries */  
     for(i=0; i < enginedata->nbnds; ++i){  
         switch(bnd_kind(enginedata->bndlist[i])){  
             case e_bnd_rel: /* real-valued boundary relation */  
                 gout[i] = bndman_real_eval(enginedata->bndlist[i]);  
 #ifdef ROOT_DEBUG  
                 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);  
                 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);  
                 ASC_FREE(relname);  
 #endif  
                 break;  
             case e_bnd_logrel:  
                 if(bndman_log_eval(enginedata->bndlist[i])){  
                     CONSOLE_DEBUG("bnd[%d] = TRUE",i);  
 #ifdef ROOT_DEBUG  
                     relname = bnd_make_name(integ->system,enginedata->bndlist[i]);  
                     CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);  
                     ASC_FREE(relname);  
 #endif  
                     gout[i] = +1.0;  
                 }else{  
                     CONSOLE_DEBUG("bnd[%d] = FALSE",i);  
                     gout[i] = -1.0;  
                 }  
                 break;  
             case e_bnd_undefined:  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");  
                 return 1;  
         }  
     }  
   
     return 0; /* no way to detect errors in bndman_*_eval at this stage */  
 }  
   
   
 /*----------------------------------------------  
   FULL JACOBIAN PRECONDITIONER -- EXPERIMENTAL.  
 */  
   
 static void integrator_ida_pcreate_jacobian(IntegratorSystem *integ){  
     IntegratorIdaData *enginedata =integ->enginedata;  
     IntegratorIdaPrecDataJacobian *precdata;  
     precdata = ASC_NEW(IntegratorIdaPrecDataJacobian);  
     mtx_matrix_t P;  
     asc_assert(integ->n_y);  
     precdata->L = linsolqr_create_default();  
   
     /* allocate matrix to be used by linsolqr */  
     P = mtx_create();  
     mtx_set_order(P, integ->n_y);  
     linsolqr_set_matrix(precdata->L, P);  
   
     enginedata->pfree = &integrator_ida_pfree_jacobian;  
     enginedata->precdata = precdata;  
     CONSOLE_DEBUG("Allocated memory for Full Jacobian preconditioner");  
 }  
   
 static void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata){  
     mtx_matrix_t P;  
     IntegratorIdaPrecDataJacobian *precdata;  
   
     if(enginedata->precdata){  
         precdata = (IntegratorIdaPrecDataJacobian *)enginedata->precdata;  
         P = linsolqr_get_matrix(precdata->L);  
         mtx_destroy(P);  
         linsolqr_destroy(precdata->L);  
         ASC_FREE(precdata);  
         enginedata->precdata = NULL;  
   
         CONSOLE_DEBUG("Freed memory for Full Jacobian preconditioner");  
     }  
     enginedata->pfree = NULL;  
 }  
   
 /**  
     EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers  
   
     'setup' function.  
 */  
 static int integrator_ida_psetup_jacobian(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          realtype c_j, void *p_data,  
          N_Vector tmp1, N_Vector tmp2,  
          N_Vector tmp3  
 ){  
     int i, j, res;  
     IntegratorSystem *integ;  
     IntegratorIdaData *enginedata;  
     IntegratorIdaPrecDataJacobian *precdata;  
     linsolqr_system_t L;  
     mtx_matrix_t P;  
     struct rel_relation **relptr;  
   
     integ = (IntegratorSystem *)p_data;  
     enginedata = integ->enginedata;  
     precdata = (IntegratorIdaPrecDataJacobian *)(enginedata->precdata);  
     double *derivatives;  
     struct var_variable **variables;  
     int count, status;  
     char *relname;  
     mtx_coord_t C;  
   
     L = precdata->L;  
     P = linsolqr_get_matrix(L);  
     mtx_clear(P);  
   
     CONSOLE_DEBUG("Setting up Jacobian preconditioner");  
   
     variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);  
     derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);  
   
     /**  
         @TODO FIXME here we are using the very inefficient and contorted approach  
         of calculating the whole jacobian, then extracting just the diagonal elements.  
     */  
   
     for(i=0, relptr = enginedata->rellist;  
             i< enginedata->nrels && relptr != NULL;  
             ++i, ++relptr  
     ){  
         /* get derivatives for this particular relation */  
         status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);  
         if(status){  
             relname = rel_make_name(integ->system, *relptr);  
             CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);  
             ASC_FREE(relname);  
             break;  
         }  
         /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */  
         /* find the diagonal elements */  
         for(j=0; j<count; ++j){  
             if(var_deriv(variables[j])){  
                 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), c_j * derivatives[j]);  
             }else{  
                 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), derivatives[j]);  
             }  
         }  
     }  
   
     mtx_assemble(P);  
   
     if(status){  
         CONSOLE_DEBUG("Error found when evaluating derivatives");  
         res = 1; goto finish; /* recoverable */  
     }  
   
     integrator_ida_write_incidence(integ);  
   
     res = 0;  
 finish:  
     ASC_FREE(variables);  
     ASC_FREE(derivatives);  
     return res;  
 };  
   
 /**  
     EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers  
   
     'solve' function.  
 */  
 static int integrator_ida_psolve_jacobian(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          N_Vector rvec, N_Vector zvec,  
          realtype c_j, realtype delta, void *p_data,  
          N_Vector tmp  
 ){  
     IntegratorSystem *integ;  
     IntegratorIdaData *data;  
     IntegratorIdaPrecDataJacobian *precdata;  
     integ = (IntegratorSystem *)p_data;  
     data = integ->enginedata;  
     precdata = (IntegratorIdaPrecDataJacobian *)(data->precdata);  
     linsolqr_system_t L = precdata->L;  
   
     linsolqr_add_rhs(L,NV_DATA_S(rvec),FALSE);  
   
     mtx_region_t R;  
     R.row.low = R.col.low = 0;  
     R.row.high = R.col.high = mtx_order(linsolqr_get_matrix(L)) - 1;  
     linsolqr_set_region(L,R);  
   
     linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));  
     linsolqr_reorder(L, &R, linsolqr_rmethod(L));  
   
     /// @TODO more here  
   
     linsolqr_remove_rhs(L,NV_DATA_S(rvec));  
   
     CONSOLE_DEBUG("Solving Jacobian preconditioner (c_j = %f)",c_j);  
     return 0;  
 };  
   
   
 /*----------------------------------------------  
   JACOBI PRECONDITIONER -- EXPERIMENTAL.  
 */  
   
 static void integrator_ida_pcreate_jacobi(IntegratorSystem *integ){  
     IntegratorIdaData *enginedata =integ->enginedata;  
     IntegratorIdaPrecDataJacobi *precdata;  
     precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);  
   
     asc_assert(integ->n_y);  
     precdata->PIii = N_VNew_Serial(integ->n_y);  
   
     enginedata->pfree = &integrator_ida_pfree_jacobi;  
     enginedata->precdata = precdata;  
     CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");  
 }  
   
 static void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){  
     if(enginedata->precdata){  
         IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;  
         N_VDestroy_Serial(precdata->PIii);  
   
         ASC_FREE(precdata);  
         enginedata->precdata = NULL;  
         CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");  
     }  
     enginedata->pfree = NULL;  
 }  
   
 /**  
     EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers  
   
     'setup' function.  
 */  
 static int integrator_ida_psetup_jacobi(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          realtype c_j, void *p_data,  
          N_Vector tmp1, N_Vector tmp2,  
          N_Vector tmp3  
 ){  
     int i, j, res;  
     IntegratorSystem *integ;  
     IntegratorIdaData *enginedata;  
     IntegratorIdaPrecDataJacobi *precdata;  
     struct rel_relation **relptr;  
   
     integ = (IntegratorSystem *)p_data;  
     enginedata = integ->enginedata;  
     precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);  
     double *derivatives;  
     struct var_variable **variables;  
     int count, status;  
     char *relname;  
   
     CONSOLE_DEBUG("Setting up Jacobi preconditioner");  
   
     variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);  
     derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);  
   
     /**  
         @TODO FIXME here we are using the very inefficient and contorted approach  
         of calculating the whole jacobian, then extracting just the diagonal elements.  
     */  
   
     for(i=0, relptr = enginedata->rellist;  
             i< enginedata->nrels && relptr != NULL;  
             ++i, ++relptr  
     ){  
   
         /* get derivatives for this particular relation */  
         status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);  
         if(status){  
             relname = rel_make_name(integ->system, *relptr);  
             CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);  
             ASC_FREE(relname);  
             break;  
         }  
         /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */  
         /* find the diagonal elements */  
         for(j=0; j<count; ++j){  
             if(var_sindex(variables[j])==i){  
                 if(var_deriv(variables[j])){  
                     NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);  
                 }else{  
                     NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];  
                 }  
   
             }  
         }  
 #ifdef PREC_DEBUG  
         CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));  
 #endif  
     }  
   
     if(status){  
         CONSOLE_DEBUG("Error found when evaluating derivatives");  
         res = 1; goto finish; /* recoverable */  
     }  
   
     integrator_ida_write_incidence(integ);  
   
     res = 0;  
 finish:  
     ASC_FREE(variables);  
     ASC_FREE(derivatives);  
     return res;  
 };  
   
 /**  
     EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers  
   
     'solve' function.  
 */  
 static int integrator_ida_psolve_jacobi(realtype tt,  
          N_Vector yy, N_Vector yp, N_Vector rr,  
          N_Vector rvec, N_Vector zvec,  
          realtype c_j, realtype delta, void *p_data,  
          N_Vector tmp  
 ){  
     IntegratorSystem *integ;  
     IntegratorIdaData *data;  
     IntegratorIdaPrecDataJacobi *precdata;  
     integ = (IntegratorSystem *)p_data;  
     data = integ->enginedata;  
     precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);  
   
     CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);  
     N_VProd(precdata->PIii, rvec, zvec);  
     return 0;  
 };  
   
1026  /*----------------------------------------------  /*----------------------------------------------
1027    STATS    STATS
1028  */  */
# Line 2245  static int integrator_ida_write_matrix(c Line 1217  static int integrator_ida_write_matrix(c
1217      This routine outputs matrix structure in a crude text format, for the sake      This routine outputs matrix structure in a crude text format, for the sake
1218      of debugging.      of debugging.
1219  */  */
1220  static void integrator_ida_write_incidence(IntegratorSystem *integ){  void integrator_ida_write_incidence(IntegratorSystem *integ){
1221      int i, j;      int i, j;
1222      struct rel_relation **relptr;      struct rel_relation **relptr;
1223      IntegratorIdaData *enginedata = integ->enginedata;      IntegratorIdaData *enginedata = integ->enginedata;
Legend:
Removed from v.2368  
changed lines
  Added in v.2369
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22