/[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 1305 by johnpye, Thu Mar 1 06:04:21 2007 UTC revision 1306 by johnpye, Sat Mar 3 11:50:47 2007 UTC
# Line 70  Line 70 
70  #include <solver/slv_client.h>  #include <solver/slv_client.h>
71  #include <solver/relman.h>  #include <solver/relman.h>
72  #include <solver/block.h>  #include <solver/block.h>
73    #include <solver/slv_stdcalls.h>
74    
75  #include "idalinear.h"  #include "idalinear.h"
76  #include "idaanalyse.h"  #include "idaanalyse.h"
# Line 114  const IntegratorInternals integrator_ida Line 115  const IntegratorInternals integrator_ida
115      ,integrator_analyse_dae /* note, this routine is back in integrator.c */      ,integrator_analyse_dae /* note, this routine is back in integrator.c */
116  #endif  #endif
117      ,integrator_ida_solve      ,integrator_ida_solve
118      ,NULL /* writematrixfn */      ,integrator_ida_write_matrix
119      ,integrator_ida_debug      ,integrator_ida_debug
120      ,integrator_ida_free      ,integrator_ida_free
121      ,INTEG_IDA      ,INTEG_IDA
# Line 129  const IntegratorInternals integrator_ida Line 130  const IntegratorInternals integrator_ida
130  struct IntegratorIdaDataStruct;  struct IntegratorIdaDataStruct;
131    
132  /* functions for allocating storage for and freeing preconditioner data */  /* functions for allocating storage for and freeing preconditioner data */
133  typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *blsys);  typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *sys);
134  typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);  typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
135    
136  /**  /**
# Line 206  static void integrator_dae_show_var(Inte Line 207  static void integrator_dae_show_var(Inte
207    
208  int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);  int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
209  void integrator_ida_write_stats(IntegratorIdaStats *stats);  void integrator_ida_write_stats(IntegratorIdaStats *stats);
210  void integrator_ida_write_incidence(IntegratorSystem *blsys);  void integrator_ida_write_incidence(IntegratorSystem *sys);
211    
212  /*------  /*------
213    Jacobi preconditioner -- experimental    Jacobi preconditioner -- experimental
# Line 226  int integrator_ida_psolve_jacobi(realtyp Line 227  int integrator_ida_psolve_jacobi(realtyp
227           N_Vector tmp           N_Vector tmp
228  );  );
229    
230  void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys);  void integrator_ida_pcreate_jacobi(IntegratorSystem *sys);
231    
232  void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);  void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);
233    
# Line 239  static const IntegratorIdaPrec prec_jaco Line 240  static const IntegratorIdaPrec prec_jaco
240  /*-------------------------------------------------------------  /*-------------------------------------------------------------
241    SETUP/TEARDOWN ROUTINES    SETUP/TEARDOWN ROUTINES
242  */  */
243  void integrator_ida_create(IntegratorSystem *blsys){  void integrator_ida_create(IntegratorSystem *sys){
244      CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");      CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
245      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
246      enginedata = ASC_NEW(IntegratorIdaData);      enginedata = ASC_NEW(IntegratorIdaData);
# Line 253  void integrator_ida_create(IntegratorSys Line 254  void integrator_ida_create(IntegratorSys
254      enginedata->rfilter.matchbits =  REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;      enginedata->rfilter.matchbits =  REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
255      enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;      enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
256    
257      blsys->enginedata = (void *)enginedata;      sys->enginedata = (void *)enginedata;
258    
259      integrator_ida_params_default(blsys);      integrator_ida_params_default(sys);
260  }  }
261    
262  void integrator_ida_free(void *enginedata){  void integrator_ida_free(void *enginedata){
# Line 279  void integrator_ida_free(void *enginedat Line 280  void integrator_ida_free(void *enginedat
280  #endif  #endif
281  }  }
282    
283  IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){  IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *sys){
284      IntegratorIdaData *d;      IntegratorIdaData *d;
285      assert(blsys!=NULL);      assert(sys!=NULL);
286      assert(blsys->enginedata!=NULL);      assert(sys->enginedata!=NULL);
287      assert(blsys->engine==INTEG_IDA);      assert(sys->engine==INTEG_IDA);
288      d = ((IntegratorIdaData *)(blsys->enginedata));      d = ((IntegratorIdaData *)(sys->enginedata));
289      return d;      return d;
290  }  }
291    
# Line 310  enum ida_parameters{ Line 311  enum ida_parameters{
311    
312  /**  /**
313      Here the full set of parameters is defined, along with upper/lower bounds,      Here the full set of parameters is defined, along with upper/lower bounds,
314      etc. The values are stuck into the blsys->params structure.      etc. The values are stuck into the sys->params structure.
315    
316      To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters      To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
317      list. Then add a slv_param_*(...) statement below to define the type, description and range      list. Then add a slv_param_*(...) statement below to define the type, description and range
# Line 318  enum ida_parameters{ Line 319  enum ida_parameters{
319    
320      @return 0 on success      @return 0 on success
321  */  */
322  int integrator_ida_params_default(IntegratorSystem *blsys){  int integrator_ida_params_default(IntegratorSystem *sys){
323      asc_assert(blsys!=NULL);      asc_assert(sys!=NULL);
324      asc_assert(blsys->engine==INTEG_IDA);      asc_assert(sys->engine==INTEG_IDA);
325      slv_parameters_t *p;      slv_parameters_t *p;
326      p = &(blsys->params);      p = &(sys->params);
327    
328      slv_destroy_parms(p);      slv_destroy_parms(p);
329    
# Line 464  typedef char *IdaFlagNameFn(int); Line 465  typedef char *IdaFlagNameFn(int);
465    
466  /* return 0 on success */  /* return 0 on success */
467  int integrator_ida_solve(  int integrator_ida_solve(
468          IntegratorSystem *blsys          IntegratorSystem *sys
469          , unsigned long start_index          , unsigned long start_index
470          , unsigned long finish_index          , unsigned long finish_index
471  ){  ){
# Line 479  int integrator_ida_solve( Line 480  int integrator_ida_solve(
480      IdaFlagNameFn *flagnamefn;      IdaFlagNameFn *flagnamefn;
481      const char *flagfntype;      const char *flagfntype;
482      char *pname = NULL;      char *pname = NULL;
483    #ifdef SOLVE_DEBUG
484      char *varname;      char *varname;
485    #endif
486      int i;      int i;
487      const IntegratorIdaPrec *prec = NULL;      const IntegratorIdaPrec *prec = NULL;
488      int icopt; /* initial conditions strategy */      int icopt; /* initial conditions strategy */
489    
490      CONSOLE_DEBUG("STARTING IDA...");      CONSOLE_DEBUG("STARTING IDA...");
491    
492      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(sys);
493    
494      enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);      enginedata->safeeval = SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_SAFEEVAL);
495      CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);      CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
496    
497      /* store reference to list of relations (in enginedata) */      /* store reference to list of relations (in enginedata) */
498      enginedata->nrels = slv_get_num_solvers_rels(blsys->system);      enginedata->nrels = slv_get_num_solvers_rels(sys->system);
499      enginedata->rellist = slv_get_solvers_rel_list(blsys->system);      enginedata->rellist = slv_get_solvers_rel_list(sys->system);
500      enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);      enginedata->bndlist = slv_get_solvers_bnd_list(sys->system);
501    
502      CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);      CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
503      CONSOLE_DEBUG("Number of dependent vars: %d",blsys->n_y);      CONSOLE_DEBUG("Number of dependent vars: %d",sys->n_y);
504      size = blsys->n_y;      size = sys->n_y;
505    
506      if(enginedata->nrels!=size){      if(enginedata->nrels!=size){
507          ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);          ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
# Line 506  int integrator_ida_solve( Line 509  int integrator_ida_solve(
509      }      }
510    
511  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
512      integrator_ida_debug(blsys,stderr);      integrator_ida_debug(sys,stderr);
513  #endif  #endif
514    
515      /* retrieve initial values from the system */      /* retrieve initial values from the system */
516    
517      /** @TODO fix this, the starting time != first sample */      /** @TODO fix this, the starting time != first sample */
518      t0 = integrator_get_t(blsys);      t0 = integrator_get_t(sys);
519      CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);      CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
520    
521      CONSOLE_DEBUG("RETRIEVING y0");      CONSOLE_DEBUG("RETRIEVING y0");
522    
523      y0 = N_VNew_Serial(size);      y0 = N_VNew_Serial(size);
524      integrator_get_y(blsys,NV_DATA_S(y0));      integrator_get_y(sys,NV_DATA_S(y0));
525    
526  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
527      CONSOLE_DEBUG("RETRIEVING yp0");      CONSOLE_DEBUG("RETRIEVING yp0");
528  #endif  #endif
529    
530      yp0 = N_VNew_Serial(size);      yp0 = N_VNew_Serial(size);
531      integrator_get_ydot(blsys,NV_DATA_S(yp0));      integrator_get_ydot(sys,NV_DATA_S(yp0));
532    
533  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
534      N_VPrint_Serial(yp0);      N_VPrint_Serial(yp0);
# Line 536  int integrator_ida_solve( Line 539  int integrator_ida_solve(
539      ida_mem = IDACreate();      ida_mem = IDACreate();
540    
541      /* relative error tolerance */      /* relative error tolerance */
542      reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);      reltol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_RTOL);
543      CONSOLE_DEBUG("rtol = %8.2e",reltol);      CONSOLE_DEBUG("rtol = %8.2e",reltol);
544    
545      /* allocate internal memory */      /* allocate internal memory */
546      if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){      if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
547          /* vector of absolute tolerances */          /* vector of absolute tolerances */
548          CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");          CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
549          abstolvect = N_VNew_Serial(size);          abstolvect = N_VNew_Serial(size);
550          integrator_get_atol(blsys,NV_DATA_S(abstolvect));          integrator_get_atol(sys,NV_DATA_S(abstolvect));
551    
552          flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);          flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
553    
554          N_VDestroy_Serial(abstolvect);          N_VDestroy_Serial(abstolvect);
555      }else{      }else{
556          /* scalar absolute tolerance (one value for all) */          /* scalar absolute tolerance (one value for all) */
557          abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);          abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
558          CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);          CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
559          flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);          flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
560      }      }
# Line 568  int integrator_ida_solve( Line 571  int integrator_ida_solve(
571      }/* else success */      }/* else success */
572    
573      /* set optional inputs... */      /* set optional inputs... */
574      IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);      IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)sys);
575      IDASetRdata(ida_mem, (void *)blsys);      IDASetRdata(ida_mem, (void *)sys);
576      IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));      IDASetMaxStep(ida_mem, integrator_get_maxstep(sys));
577      IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));      IDASetInitStep(ida_mem, integrator_get_stepzero(sys));
578      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(sys));
579      if(integrator_get_minstep(blsys)>0){      if(integrator_get_minstep(sys)>0){
580          ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");          ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
581      }      }
582    
583      CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));      CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXNCF));
584      IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));      IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXNCF));
585    
586      CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXORD));      CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXORD));
587      IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXORD));      IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXORD));
588    
589      /* there's no capability for setting *minimum* step size in IDA */      /* there's no capability for setting *minimum* step size in IDA */
590    
591    
592      /* attach linear solver module, using the default value of maxl */      /* attach linear solver module, using the default value of maxl */
593      linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);      linsolver = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_LINSOLVER);
594      CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);      CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
595      if(strcmp(linsolver,"ASCEND")==0){      if(strcmp(linsolver,"ASCEND")==0){
596          CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",size);          CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",size);
597          IDAASCEND(ida_mem,size);          IDAASCEND(ida_mem,size);
598          IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)blsys);          IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)sys);
599    
600          flagfntype = "IDAASCEND";          flagfntype = "IDAASCEND";
601          flagfn = &IDAASCENDGetLastFlag;          flagfn = &IDAASCENDGetLastFlag;
# Line 609  int integrator_ida_solve( Line 612  int integrator_ida_solve(
612              default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;              default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
613          }          }
614    
615          if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){          if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
616              CONSOLE_DEBUG("USING AUTODIFF");              CONSOLE_DEBUG("USING AUTODIFF");
617              flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);              flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)sys);
618              switch(flag){              switch(flag){
619                  case IDADENSE_SUCCESS: break;                  case IDADENSE_SUCCESS: break;
620                  default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;                  default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
# Line 627  int integrator_ida_solve( Line 630  int integrator_ida_solve(
630          /* remaining methods are all SPILS */          /* remaining methods are all SPILS */
631          CONSOLE_DEBUG("IDA SPILS");          CONSOLE_DEBUG("IDA SPILS");
632    
633          maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);          maxl = SLV_PARAM_INT(&(sys->params),IDA_PARAM_MAXL);
634          CONSOLE_DEBUG("maxl = %d",maxl);          CONSOLE_DEBUG("maxl = %d",maxl);
635    
636          /* what preconditioner? */          /* what preconditioner? */
637          pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);          pname = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_PREC);
638          if(strcmp(pname,"NONE")==0){          if(strcmp(pname,"NONE")==0){
639              prec = NULL;              prec = NULL;
640          }else if(strcmp(pname,"JACOBI")==0){          }else if(strcmp(pname,"JACOBI")==0){
# Line 658  int integrator_ida_solve( Line 661  int integrator_ida_solve(
661    
662          if(prec){          if(prec){
663              /* assign the preconditioner to the linear solver */              /* assign the preconditioner to the linear solver */
664              (prec->pcreate)(blsys);              (prec->pcreate)(sys);
665              IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)blsys);              IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)sys);
666              CONSOLE_DEBUG("PRECONDITIONER = %s",pname);              CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
667          }else{          }else{
668              CONSOLE_DEBUG("No preconditioner");              CONSOLE_DEBUG("No preconditioner");
# Line 678  int integrator_ida_solve( Line 681  int integrator_ida_solve(
681          }/* else success */          }/* else success */
682    
683          /* assign the J*v function */          /* assign the J*v function */
684          if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){          if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
685              CONSOLE_DEBUG("USING AUTODIFF");              CONSOLE_DEBUG("USING AUTODIFF");
686              flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);              flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)sys);
687              if(flag==IDASPILS_MEM_NULL){              if(flag==IDASPILS_MEM_NULL){
688                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
689                  return 10;                  return 10;
# Line 694  int integrator_ida_solve( Line 697  int integrator_ida_solve(
697    
698          if(strcmp(linsolver,"SPGMR")==0){          if(strcmp(linsolver,"SPGMR")==0){
699              /* select Gram-Schmidt orthogonalisation */              /* select Gram-Schmidt orthogonalisation */
700              if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){              if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_GSMODIFIED)){
701                  CONSOLE_DEBUG("USING MODIFIED GS");                  CONSOLE_DEBUG("USING MODIFIED GS");
702                  flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);                  flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
703                  if(flag!=IDASPILS_SUCCESS){                  if(flag!=IDASPILS_SUCCESS){
# Line 718  int integrator_ida_solve( Line 721  int integrator_ida_solve(
721    
722      /* calculate initial conditions */      /* calculate initial conditions */
723      icopt = 0;      icopt = 0;
724      if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"Y")==0){      if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"Y")==0){
725          CONSOLE_DEBUG("Solving initial conditions using values of yddot");          CONSOLE_DEBUG("Solving initial conditions using values of yddot");
726          icopt = IDA_Y_INIT;          icopt = IDA_Y_INIT;
727          asc_assert(icopt!=0);          asc_assert(icopt!=0);
728      }else if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){      }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
729          CONSOLE_DEBUG("Solving initial conditions using values of yd");          CONSOLE_DEBUG("Solving initial conditions using values of yd");
730          icopt = IDA_YA_YDP_INIT;          icopt = IDA_YA_YDP_INIT;
731          asc_assert(icopt!=0);          asc_assert(icopt!=0);
732          id = N_VNew_Serial(blsys->n_y);          id = N_VNew_Serial(sys->n_y);
733          for(i=0; i < blsys->n_y; ++i){          for(i=0; i < sys->n_y; ++i){
734              if(blsys->ydot[i] == NULL){              if(sys->ydot[i] == NULL){
735                  NV_Ith_S(id,i) = 0.0;                  NV_Ith_S(id,i) = 0.0;
736  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
737                  varname = var_make_name(blsys->system,blsys->y[i]);                  varname = var_make_name(sys->system,sys->y[i]);
738                  CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);                  CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
739                  ASC_FREE(varname);                  ASC_FREE(varname);
740  #endif  #endif
# Line 744  int integrator_ida_solve( Line 747  int integrator_ida_solve(
747          }          }
748          IDASetId(ida_mem, id);          IDASetId(ida_mem, id);
749          N_VDestroy_Serial(id);          N_VDestroy_Serial(id);
750      }else if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"NONE")==0){      }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"NONE")==0){
751          ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");          ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
752      }else{      }else{
753          ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");          ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
754      }      }
755    
756      if(icopt){      if(icopt){
757          blsys->currentstep=0;          sys->currentstep=0;
758          t_index=start_index + 1;          t_index=start_index + 1;
759          tout1 = samplelist_get(blsys->samples, t_index);          tout1 = samplelist_get(sys->samples, t_index);
760    
761          CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);          CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
762    
# Line 824  int integrator_ida_solve( Line 827  int integrator_ida_solve(
827      /* optionally, specify ROO-FINDING PROBLEM */      /* optionally, specify ROO-FINDING PROBLEM */
828    
829      /* -- set up the IntegratorReporter */      /* -- set up the IntegratorReporter */
830      integrator_output_init(blsys);      integrator_output_init(sys);
831    
832      /* -- store the initial values of all the stuff */      /* -- store the initial values of all the stuff */
833      integrator_output_write(blsys);      integrator_output_write(sys);
834      integrator_output_write_obs(blsys);      integrator_output_write_obs(sys);
835    
836      /* specify where the returned values should be stored */      /* specify where the returned values should be stored */
837      yret = y0;      yret = y0;
838      ypret = yp0;      ypret = yp0;
839    
840      /* advance solution in time, return values as yret and derivatives as ypret */      /* advance solution in time, return values as yret and derivatives as ypret */
841      blsys->currentstep=1;      sys->currentstep=1;
842      for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){      for(t_index=start_index+1;t_index <= finish_index;++t_index, ++sys->currentstep){
843          t = samplelist_get(blsys->samples, t_index);          t = samplelist_get(sys->samples, t_index);
844          t0 = integrator_get_t(blsys);          t0 = integrator_get_t(sys);
845          asc_assert(t > t0);          asc_assert(t > t0);
846    
847  #ifdef SOLVE_DEBUG  #ifdef SOLVE_DEBUG
# Line 848  int integrator_ida_solve( Line 851  int integrator_ida_solve(
851          flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);          flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
852    
853          /* pass the values of everything back to the compiler */          /* pass the values of everything back to the compiler */
854          integrator_set_t(blsys, (double)tret);          integrator_set_t(sys, (double)tret);
855          integrator_set_y(blsys, NV_DATA_S(yret));          integrator_set_y(sys, NV_DATA_S(yret));
856          integrator_set_ydot(blsys, NV_DATA_S(ypret));          integrator_set_ydot(sys, NV_DATA_S(ypret));
857    
858          if(flag<0){          if(flag<0){
859              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
860              break;              break;
861          }          }
862    
863          /* -- do something so that blsys knows the values of tret, yret and ypret */          /* -- do something so that sys knows the values of tret, yret and ypret */
864    
865          /* -- store the current values of all the stuff */          /* -- store the current values of all the stuff */
866          integrator_output_write(blsys);          integrator_output_write(sys);
867          integrator_output_write_obs(blsys);          integrator_output_write_obs(sys);
868      }      }
869    
870      /* -- close the IntegratorReporter */      /* -- close the IntegratorReporter */
871      integrator_output_close(blsys);      integrator_output_close(sys);
872    
873      /* get optional outputs */      /* get optional outputs */
874  #ifdef STATS_DEBUG  #ifdef STATS_DEBUG
# Line 936  void integrator_ida_sig(int sig){ Line 939  void integrator_ida_sig(int sig){
939      @param yy current values of dependent variable vector      @param yy current values of dependent variable vector
940      @param yp current values of derivatives of dependent variables      @param yp current values of derivatives of dependent variables
941      @param rr the output residual vector (is we're returning data to)      @param rr the output residual vector (is we're returning data to)
942      @param res_data pointer to our stuff (blsys in this case).      @param res_data pointer to our stuff (sys in this case).
943    
944      @return 0 on success, positive on recoverable error, and      @return 0 on success, positive on recoverable error, and
945          negative on unrecoverable error.          negative on unrecoverable error.
946  */  */
947  int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){  int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
948      IntegratorSystem *blsys;      IntegratorSystem *sys;
949      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
950      int i, calc_ok, is_error;      int i, calc_ok, is_error;
951      struct rel_relation** relptr;      struct rel_relation** relptr;
# Line 953  int integrator_ida_fex(realtype tt, N_Ve Line 956  int integrator_ida_fex(realtype tt, N_Ve
956      char diffname[30];      char diffname[30];
957  #endif  #endif
958    
959      blsys = (IntegratorSystem *)res_data;      sys = (IntegratorSystem *)res_data;
960      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(sys);
961    
962  #ifdef FEX_DEBUG  #ifdef FEX_DEBUG
963      /* fprintf(stderr,"\n\n"); */      /* fprintf(stderr,"\n\n"); */
# Line 967  int integrator_ida_fex(realtype tt, N_Ve Line 970  int integrator_ida_fex(realtype tt, N_Ve
970      }      }
971    
972      /* pass the values of everything back to the compiler */      /* pass the values of everything back to the compiler */
973      integrator_set_t(blsys, (double)tt);      integrator_set_t(sys, (double)tt);
974      integrator_set_y(blsys, NV_DATA_S(yy));      integrator_set_y(sys, NV_DATA_S(yy));
975      integrator_set_ydot(blsys, NV_DATA_S(yp));      integrator_set_ydot(sys, NV_DATA_S(yp));
976    
977      /* perform bounds checking on all variables */      /* perform bounds checking on all variables */
978      if(slv_check_bounds(blsys->system, 0, -1)){      if(slv_check_bounds(sys->system, 0, -1, NULL)){
979          /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */          /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
980          return 1;          return 1;
981      }      }
# Line 1011  int integrator_ida_fex(realtype tt, N_Ve Line 1014  int integrator_ida_fex(realtype tt, N_Ve
1014    
1015              NV_Ith_S(rr,i) = resid;              NV_Ith_S(rr,i) = resid;
1016              if(!calc_ok){              if(!calc_ok){
1017                  relname = rel_make_name(blsys->system, *relptr);                  relname = rel_make_name(sys->system, *relptr);
1018                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1019                  ASC_FREE(relname);                  ASC_FREE(relname);
1020                  /* presumable some output already made? */                  /* presumable some output already made? */
# Line 1050  int integrator_ida_fex(realtype tt, N_Ve Line 1053  int integrator_ida_fex(realtype tt, N_Ve
1053    
1054  #ifdef ASC_SIGNAL_TRAPS  #ifdef ASC_SIGNAL_TRAPS
1055      }else{      }else{
1056          relname = rel_make_name(blsys->system, *relptr);          relname = rel_make_name(sys->system, *relptr);
1057          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1058          ASC_FREE(relname);          ASC_FREE(relname);
1059          is_error = 1;          is_error = 1;
# Line 1068  int integrator_ida_fex(realtype tt, N_Ve Line 1071  int integrator_ida_fex(realtype tt, N_Ve
1071      /* output residuals to console */      /* output residuals to console */
1072      CONSOLE_DEBUG("RESIDUAL OUTPUT");      CONSOLE_DEBUG("RESIDUAL OUTPUT");
1073      fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");      fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1074      for(i=0; i<blsys->n_y; ++i){      for(i=0; i<sys->n_y; ++i){
1075          varname = var_make_name(blsys->system,blsys->y[i]);          varname = var_make_name(sys->system,sys->y[i]);
1076          fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));          fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1077          if(blsys->ydot[i]){          if(sys->ydot[i]){
1078              varname = var_make_name(blsys->system,blsys->ydot[i]);              varname = var_make_name(sys->system,sys->ydot[i]);
1079              fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));              fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1080          }else{          }else{
1081              snprintf(diffname,99,"diff(%s)",varname);              snprintf(diffname,99,"diff(%s)",varname);
1082              fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));              fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1083          }          }
1084          ASC_FREE(varname);          ASC_FREE(varname);
1085          relname = rel_make_name(blsys->system,enginedata->rellist[i]);          relname = rel_make_name(sys->system,enginedata->rellist[i]);
1086          fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);          fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1087      }      }
1088  #endif  #endif
# Line 1103  int integrator_ida_djex(long int Neq, re Line 1106  int integrator_ida_djex(long int Neq, re
1106          , realtype c_j, void *jac_data, DenseMat Jac          , realtype c_j, void *jac_data, DenseMat Jac
1107          , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3          , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1108  ){  ){
1109      IntegratorSystem *blsys;      IntegratorSystem *sys;
1110      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
1111      char *relname;      char *relname;
1112  #ifdef DJEX_DEBUG  #ifdef DJEX_DEBUG
# Line 1117  int integrator_ida_djex(long int Neq, re Line 1120  int integrator_ida_djex(long int Neq, re
1120      int count, j;      int count, j;
1121      int status, is_error = 0;      int status, is_error = 0;
1122    
1123      blsys = (IntegratorSystem *)jac_data;      sys = (IntegratorSystem *)jac_data;
1124      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(sys);
1125    
1126      /* allocate space for returns from relman_diff3 */      /* allocate space for returns from relman_diff3 */
1127      /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */      /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
# Line 1126  int integrator_ida_djex(long int Neq, re Line 1129  int integrator_ida_djex(long int Neq, re
1129      derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);      derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1130    
1131      /* pass the values of everything back to the compiler */      /* pass the values of everything back to the compiler */
1132      integrator_set_t(blsys, (double)tt);      integrator_set_t(sys, (double)tt);
1133      integrator_set_y(blsys, NV_DATA_S(yy));      integrator_set_y(sys, NV_DATA_S(yy));
1134      integrator_set_ydot(blsys, NV_DATA_S(yp));      integrator_set_ydot(sys, NV_DATA_S(yp));
1135    
1136      /* perform bounds checking on all variables */      /* perform bounds checking on all variables */
1137      if(slv_check_bounds(blsys->system, 0, -1)){      if(slv_check_bounds(sys->system, 0, -1, NULL)){
1138          /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */          /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1139          return 1;          return 1;
1140      }      }
1141    
1142  #ifdef DJEX_DEBUG  #ifdef DJEX_DEBUG
1143      varlist = slv_get_solvers_var_list(blsys->system);      varlist = slv_get_solvers_var_list(sys->system);
1144    
1145      /* print vars */      /* print vars */
1146      for(i=0; i < blsys->n_y; ++i){      for(i=0; i < sys->n_y; ++i){
1147          varname = var_make_name(blsys->system, blsys->y[i]);          varname = var_make_name(sys->system, sys->y[i]);
1148          CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));          CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1149          asc_assert(NV_Ith_S(yy,i) == var_value(blsys->y[i]));          asc_assert(NV_Ith_S(yy,i) == var_value(sys->y[i]));
1150          ASC_FREE(varname);          ASC_FREE(varname);
1151      }      }
1152    
1153      /* print derivatives */      /* print derivatives */
1154      for(i=0; i < blsys->n_y; ++i){      for(i=0; i < sys->n_y; ++i){
1155          if(blsys->ydot[i]){          if(sys->ydot[i]){
1156              varname = var_make_name(blsys->system, blsys->ydot[i]);              varname = var_make_name(sys->system, sys->ydot[i]);
1157              CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));              CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(sys->ydot[i]));
1158              ASC_FREE(varname);              ASC_FREE(varname);
1159          }else{          }else{
1160              varname = var_make_name(blsys->system, blsys->y[i]);              varname = var_make_name(sys->system, sys->y[i]);
1161              CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));              CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1162              ASC_FREE(varname);              ASC_FREE(varname);
1163          }          }
# Line 1174  int integrator_ida_djex(long int Neq, re Line 1177  int integrator_ida_djex(long int Neq, re
1177          status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);          status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1178    
1179          if(status){          if(status){
1180              relname = rel_make_name(blsys->system, *relptr);              relname = rel_make_name(sys->system, *relptr);
1181              CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);              CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1182              ASC_FREE(relname);              ASC_FREE(relname);
1183              is_error = 1;              is_error = 1;
# Line 1183  int integrator_ida_djex(long int Neq, re Line 1186  int integrator_ida_djex(long int Neq, re
1186    
1187          /* output what's going on here ... */          /* output what's going on here ... */
1188  #ifdef DJEX_DEBUG  #ifdef DJEX_DEBUG
1189          relname = rel_make_name(blsys->system, *relptr);          relname = rel_make_name(sys->system, *relptr);
1190          fprintf(stderr,"%d: '%s': ",i,relname);          fprintf(stderr,"%d: '%s': ",i,relname);
1191          for(j=0;j<count;++j){          for(j=0;j<count;++j){
1192              varname = var_make_name(blsys->system, variables[j]);              varname = var_make_name(sys->system, variables[j]);
1193              if(var_deriv(variables[j])){              if(var_deriv(variables[j])){
1194                  fprintf(stderr,"  '%s'=",varname);                  fprintf(stderr,"  '%s'=",varname);
1195                  fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(blsys,variables[j]));                  fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(sys,variables[j]));
1196              }else{              }else{
1197                  fprintf(stderr,"  '%s'=y[%d]",varname,var_sindex(variables[j]));                  fprintf(stderr,"  '%s'=y[%d]",varname,var_sindex(variables[j]));
1198              }              }
# Line 1202  int integrator_ida_djex(long int Neq, re Line 1205  int integrator_ida_djex(long int Neq, re
1205          /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */          /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1206          for(j=0; j < count; ++j){          for(j=0; j < count; ++j){
1207  #ifdef DJEX_DEBUG  #ifdef DJEX_DEBUG
1208              varname = var_make_name(blsys->system,variables[j]);              varname = var_make_name(sys->system,variables[j]);
1209              fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);              fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1210              ASC_FREE(varname);              ASC_FREE(varname);
1211  #endif  #endif
# Line 1214  int integrator_ida_djex(long int Neq, re Line 1217  int integrator_ida_djex(long int Neq, re
1217  #endif  #endif
1218                  DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];                  DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1219              }else{              }else{
1220                  DENSE_ELEM(Jac,i,integrator_ida_diffindex(blsys,variables[j])) += derivatives[j] * c_j;                  DENSE_ELEM(Jac,i,integrator_ida_diffindex(sys,variables[j])) += derivatives[j] * c_j;
1221  #ifdef DJEX_DEBUG  #ifdef DJEX_DEBUG
1222                  fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);                  fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1223  #endif  #endif
# Line 1226  int integrator_ida_djex(long int Neq, re Line 1229  int integrator_ida_djex(long int Neq, re
1229      ASC_FREE(relname);      ASC_FREE(relname);
1230      CONSOLE_DEBUG("PRINTING JAC");      CONSOLE_DEBUG("PRINTING JAC");
1231      fprintf(stderr,"\t");      fprintf(stderr,"\t");
1232      for(j=0; j < blsys->n_y; ++j){      for(j=0; j < sys->n_y; ++j){
1233          if(j)fprintf(stderr,"\t");          if(j)fprintf(stderr,"\t");
1234          varname = var_make_name(blsys->system,blsys->y[j]);          varname = var_make_name(sys->system,sys->y[j]);
1235          fprintf(stderr,"%11s",varname);          fprintf(stderr,"%11s",varname);
1236          ASC_FREE(varname);          ASC_FREE(varname);
1237      }      }
1238      fprintf(stderr,"\n");      fprintf(stderr,"\n");
1239      for(i=0; i < enginedata->nrels; ++i){      for(i=0; i < enginedata->nrels; ++i){
1240          relname = rel_make_name(blsys->system, enginedata->rellist[i]);          relname = rel_make_name(sys->system, enginedata->rellist[i]);
1241          fprintf(stderr,"%s\t",relname);          fprintf(stderr,"%s\t",relname);
1242          ASC_FREE(relname);          ASC_FREE(relname);
1243    
1244          for(j=0; j < blsys->n_y; ++j){          for(j=0; j < sys->n_y; ++j){
1245              if(j)fprintf(stderr,"\t");              if(j)fprintf(stderr,"\t");
1246              fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));              fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1247          }          }
# Line 1249  int integrator_ida_djex(long int Neq, re Line 1252  int integrator_ida_djex(long int Neq, re
1252      /* test for NANs */      /* test for NANs */
1253      if(!is_error){      if(!is_error){
1254          for(i=0;i< enginedata->nrels; ++i){          for(i=0;i< enginedata->nrels; ++i){
1255              for(j=0;j<blsys->n_y;++j){              for(j=0;j<sys->n_y;++j){
1256                  if(isnan(DENSE_ELEM(Jac,i,j))){                  if(isnan(DENSE_ELEM(Jac,i,j))){
1257                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1258                      is_error=1;                      is_error=1;
# Line 1263  int integrator_ida_djex(long int Neq, re Line 1266  int integrator_ida_djex(long int Neq, re
1266  #endif  #endif
1267      }      }
1268    
1269  /*  if(integrator_ida_check_diffindex(blsys)){  /*  if(integrator_ida_check_diffindex(sys)){
1270          is_error = 1;          is_error = 1;
1271      }*/      }*/
1272    
# Line 1291  int integrator_ida_djex(long int Neq, re Line 1294  int integrator_ida_djex(long int Neq, re
1294      @param v  the vector by which the Jacobian must be multiplied to the right.      @param v  the vector by which the Jacobian must be multiplied to the right.
1295      @param Jv the output vector computed      @param Jv the output vector computed
1296      @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).      @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1297      @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)      @param jac_data pointer to our stuff (sys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1298      @param tmp1 @see tmp2      @param tmp1 @see tmp2
1299      @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.      @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.
1300      @return 0 on success      @return 0 on success
# Line 1300  int integrator_ida_jvex(realtype tt, N_V Line 1303  int integrator_ida_jvex(realtype tt, N_V
1303          , N_Vector v, N_Vector Jv, realtype c_j          , N_Vector v, N_Vector Jv, realtype c_j
1304          , void *jac_data, N_Vector tmp1, N_Vector tmp2          , void *jac_data, N_Vector tmp1, N_Vector tmp2
1305  ){  ){
1306      IntegratorSystem *blsys;      IntegratorSystem *sys;
1307      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
1308      int i, j, is_error=0;      int i, j, is_error=0;
1309      struct rel_relation** relptr;      struct rel_relation** relptr;
# Line 1317  int integrator_ida_jvex(realtype tt, N_V Line 1320  int integrator_ida_jvex(realtype tt, N_V
1320      CONSOLE_DEBUG("EVALUATING JACOBIAN...");      CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1321  #endif  #endif
1322    
1323      blsys = (IntegratorSystem *)jac_data;      sys = (IntegratorSystem *)jac_data;
1324      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(sys);
1325      varlist = slv_get_solvers_var_list(blsys->system);      varlist = slv_get_solvers_var_list(sys->system);
1326    
1327      /* pass the values of everything back to the compiler */      /* pass the values of everything back to the compiler */
1328      integrator_set_t(blsys, (double)tt);      integrator_set_t(sys, (double)tt);
1329      integrator_set_y(blsys, NV_DATA_S(yy));      integrator_set_y(sys, NV_DATA_S(yy));
1330      integrator_set_ydot(blsys, NV_DATA_S(yp));      integrator_set_ydot(sys, NV_DATA_S(yp));
1331      /* no real use for residuals (rr) here, I don't think? */      /* no real use for residuals (rr) here, I don't think? */
1332    
1333      /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */      /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
# Line 1354  int integrator_ida_jvex(realtype tt, N_V Line 1357  int integrator_ida_jvex(realtype tt, N_V
1357  #endif  #endif
1358    
1359              if(status){              if(status){
1360                  relname = rel_make_name(blsys->system, *relptr);                  relname = rel_make_name(sys->system, *relptr);
1361                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1362                  ASC_FREE(relname);                  ASC_FREE(relname);
1363                  is_error = 1;                  is_error = 1;
# Line 1369  int integrator_ida_jvex(realtype tt, N_V Line 1372  int integrator_ida_jvex(realtype tt, N_V
1372    
1373              Jv_i = 0;              Jv_i = 0;
1374              for(j=0; j < count; ++j){              for(j=0; j < count; ++j){
1375                  /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);                  /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], sys->n_y);
1376                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);                  varname = var_make_name(sys->system, enginedata->varlist[variables[j]]);
1377                  if(varname){                  if(varname){
1378                      CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);                      CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1379                      ASC_FREE(varname);                      ASC_FREE(varname);
# Line 1381  int integrator_ida_jvex(realtype tt, N_V Line 1384  int integrator_ida_jvex(realtype tt, N_V
1384    
1385                  /* we don't calculate derivatives wrt indep var */                  /* we don't calculate derivatives wrt indep var */
1386                  asc_assert(variables[j]>=0);                  asc_assert(variables[j]>=0);
1387                  if(variables[j] == blsys->x) continue;                  if(variables[j] == sys->x) continue;
1388  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
1389                  CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));                  CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1390  #endif  #endif
1391                  if(var_deriv(variables[j])){                  if(var_deriv(variables[j])){
1392  #define DIFFINDEX integrator_ida_diffindex(blsys,variables[j])  #define DIFFINDEX integrator_ida_diffindex(sys,variables[j])
1393  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
1394                      fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i                      fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1395                          , derivatives[j] * NV_Ith_S(v,DIFFINDEX)                          , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
# Line 1394  int integrator_ida_jvex(realtype tt, N_V Line 1397  int integrator_ida_jvex(realtype tt, N_V
1397                          , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)                          , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1398                      );                      );
1399  #endif  #endif
1400                      asc_assert(blsys->ydot[DIFFINDEX]==variables[j]);                      asc_assert(sys->ydot[DIFFINDEX]==variables[j]);
1401                      Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;                      Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1402  #undef DIFFINDEX  #undef DIFFINDEX
1403                  }else{                  }else{
1404  #define VARINDEX var_sindex(variables[j])  #define VARINDEX var_sindex(variables[j])
1405  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
1406                      asc_assert(blsys->y[VARINDEX]==variables[j]);                      asc_assert(sys->y[VARINDEX]==variables[j]);
1407                      fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"                      fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1408                          , i                          , i
1409                          , derivatives[j] * NV_Ith_S(v,VARINDEX)                          , derivatives[j] * NV_Ith_S(v,VARINDEX)
# Line 1416  int integrator_ida_jvex(realtype tt, N_V Line 1419  int integrator_ida_jvex(realtype tt, N_V
1419              NV_Ith_S(Jv,i) = Jv_i;              NV_Ith_S(Jv,i) = Jv_i;
1420  #ifdef JEX_DEBUG  #ifdef JEX_DEBUG
1421              CONSOLE_DEBUG("rel = %p",*relptr);              CONSOLE_DEBUG("rel = %p",*relptr);
1422              relname = rel_make_name(blsys->system, *relptr);              relname = rel_make_name(sys->system, *relptr);
1423              CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));              CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1424              ASC_FREE(relname);              ASC_FREE(relname);
1425              return 1;              return 1;
# Line 1424  int integrator_ida_jvex(realtype tt, N_V Line 1427  int integrator_ida_jvex(realtype tt, N_V
1427          }          }
1428  #ifdef ASC_SIGNAL_TRAPS  #ifdef ASC_SIGNAL_TRAPS
1429      }else{      }else{
1430          relname = rel_make_name(blsys->system, *relptr);          relname = rel_make_name(sys->system, *relptr);
1431          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1432          ASC_FREE(relname);          ASC_FREE(relname);
1433          is_error = 1;          is_error = 1;
# Line 1453  int integrator_ida_sjex(long int Neq, re Line 1456  int integrator_ida_sjex(long int Neq, re
1456    JACOBI PRECONDITIONER -- EXPERIMENTAL.    JACOBI PRECONDITIONER -- EXPERIMENTAL.
1457  */  */
1458    
1459  void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){  void integrator_ida_pcreate_jacobi(IntegratorSystem *sys){
1460      IntegratorIdaData * enginedata =blsys->enginedata;      IntegratorIdaData * enginedata =sys->enginedata;
1461      IntegratorIdaPrecDataJacobi *precdata;      IntegratorIdaPrecDataJacobi *precdata;
1462      precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);      precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1463    
1464      asc_assert(blsys->n_y);      asc_assert(sys->n_y);
1465      precdata->PIii = N_VNew_Serial(blsys->n_y);      precdata->PIii = N_VNew_Serial(sys->n_y);
1466    
1467      enginedata->pfree = &integrator_ida_pfree_jacobi;      enginedata->pfree = &integrator_ida_pfree_jacobi;
1468      enginedata->precdata = precdata;      enginedata->precdata = precdata;
# Line 1490  int integrator_ida_psetup_jacobi(realtyp Line 1493  int integrator_ida_psetup_jacobi(realtyp
1493           N_Vector tmp3           N_Vector tmp3
1494  ){  ){
1495      int i, j, res;      int i, j, res;
1496      IntegratorSystem *blsys;      IntegratorSystem *sys;
1497      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
1498      IntegratorIdaPrecDataJacobi *precdata;      IntegratorIdaPrecDataJacobi *precdata;
1499      struct rel_relation **relptr;      struct rel_relation **relptr;
1500    
1501      blsys = (IntegratorSystem *)p_data;      sys = (IntegratorSystem *)p_data;
1502      enginedata = blsys->enginedata;      enginedata = sys->enginedata;
1503      precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);      precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1504      double *derivatives;      double *derivatives;
1505      struct var_variable **variables;      struct var_variable **variables;
# Line 1521  int integrator_ida_psetup_jacobi(realtyp Line 1524  int integrator_ida_psetup_jacobi(realtyp
1524          /* get derivatives for this particular relation */          /* get derivatives for this particular relation */
1525          status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);          status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1526          if(status){          if(status){
1527              relname = rel_make_name(blsys->system, *relptr);              relname = rel_make_name(sys->system, *relptr);
1528              CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);              CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1529              ASC_FREE(relname);              ASC_FREE(relname);
1530              break;              break;
# Line 1548  int integrator_ida_psetup_jacobi(realtyp Line 1551  int integrator_ida_psetup_jacobi(realtyp
1551          res = 1; goto finish; /* recoverable */          res = 1; goto finish; /* recoverable */
1552      }      }
1553    
1554      integrator_ida_write_incidence(blsys);      integrator_ida_write_incidence(sys);
1555    
1556      res = 0;      res = 0;
1557  finish:  finish:
# Line 1568  int integrator_ida_psolve_jacobi(realtyp Line 1571  int integrator_ida_psolve_jacobi(realtyp
1571           realtype c_j, realtype delta, void *p_data,           realtype c_j, realtype delta, void *p_data,
1572           N_Vector tmp           N_Vector tmp
1573  ){  ){
1574      IntegratorSystem *blsys;      IntegratorSystem *sys;
1575      IntegratorIdaData *data;      IntegratorIdaData *data;
1576      IntegratorIdaPrecDataJacobi *precdata;      IntegratorIdaPrecDataJacobi *precdata;
1577      blsys = (IntegratorSystem *)p_data;      sys = (IntegratorSystem *)p_data;
1578      data = blsys->enginedata;      data = sys->enginedata;
1579      precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);      precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1580    
1581      CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);      CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
# Line 1624  enum integrator_ida_write_jac_enum{ Line 1627  enum integrator_ida_write_jac_enum{
1627  };  };
1628    
1629  /**  /**
1630      @TODO COMPLETE THIS...      Our tasks here is to write the matrices that IDA *should* be seeing. We
1631        are actually making calls to relman_diff in order to do that, so we're
1632        really going back to the variables in the actualy system and computing
1633        row by row what the values are. This should mean just a single call to
1634        each blackbox present in the system (if blackbox caching is working
1635        correctly).
1636  */  */
1637  void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){  int integrator_ida_write_matrix(const IntegratorSystem *sys, FILE *f, const char *type){
1638      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
1639      MM_typecode matcode;      MM_typecode matcode;
1640      int nnz, rhomax;      int nnz, rhomax;
# Line 1634  void integrator_ida_write_jacobian(Integ Line 1642  void integrator_ida_write_jacobian(Integ
1642      struct var_variable **variables;      struct var_variable **variables;
1643      struct rel_relation **relptr;      struct rel_relation **relptr;
1644      int i, j, status, count;      int i, j, status, count;
1645      char *relname;      char *relname, *varname;
1646        enum integrator_ida_write_jac_enum type1;
1647      var_filter_t vfilter = {      var_filter_t vfilter = {
1648          VAR_SVAR | VAR_FIXED | VAR_DERIV          VAR_SVAR | VAR_FIXED | VAR_DERIV
1649         ,VAR_SVAR | 0         | 0         ,VAR_SVAR | 0         | 0
1650      };      };
1651      enginedata = (IntegratorIdaData *)blsys->enginedata;  
1652        if(NULL==type || 0 != strcmp(type,"ydot")){
1653            type1 = II_WRITE_Y;
1654            vfilter.matchvalue &= ~VAR_DERIV; /* clear the VAR_DERIV bit */
1655        }else{
1656            type1 = II_WRITE_YDOT;
1657            vfilter.matchvalue |= VAR_DERIV; /* set the VAR_DERIV bit */
1658        }
1659    
1660        enginedata = (IntegratorIdaData *)sys->enginedata;
1661        if(!enginedata->rellist || !enginedata->nrels){
1662            enginedata->nrels = slv_get_num_solvers_rels(sys->system);
1663            enginedata->rellist = slv_get_solvers_rel_list(sys->system);
1664        }
1665    
1666      /* number of non-zeros for all the non-FIXED solver_vars,      /* number of non-zeros for all the non-FIXED solver_vars,
1667          in all the active included equality relations.          in all the active included equality relations.
1668      */      */
1669      nnz = relman_jacobian_count(enginedata->rellist, enginedata->nrels      nnz = relman_jacobian_count(enginedata->rellist, sys->n_y
1670          , &enginedata->vfilter, &enginedata->rfilter          , &vfilter, &enginedata->rfilter
1671          , &rhomax          , &rhomax /* = the bandwidth */
1672      );      );
1673    
1674      /* we must have found the same number of relations */      CONSOLE_DEBUG("Found %d non-zeros",nnz);
     asc_assert(rhomax == enginedata->nrels);  
1675    
1676      /* output the mmio file header, now that we know our size*/      /* output the mmio file header, now that we know our size*/
1677      /* note that we are asserting that our problem is square */      /* note that we are asserting that our problem is square */
# Line 1662  void integrator_ida_write_jacobian(Integ Line 1682  void integrator_ida_write_jacobian(Integ
1682      mm_write_banner(f, matcode);      mm_write_banner(f, matcode);
1683      mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);      mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1684    
1685      variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);      variables = ASC_NEW_ARRAY(struct var_variable *, sys->n_y * 2);
1686      derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);      derivatives = ASC_NEW_ARRAY(double, sys->n_y * 2);
1687    
1688      CONSOLE_DEBUG("Writing sparse Jacobian to file...");      CONSOLE_DEBUG("Writing matrix dF/d%s (%d rels) to file...",(type1==II_WRITE_Y ? "y" : "y'"),enginedata->nrels);
1689    
1690      for(i=0, relptr = enginedata->rellist;      for(i=0, relptr = enginedata->rellist;
1691              i< enginedata->nrels && relptr != NULL;              i< enginedata->nrels && relptr != NULL;
1692              ++i, ++relptr              ++i, ++relptr
1693      ){      ){
1694          relname = rel_make_name(blsys->system, *relptr);          relname = rel_make_name(sys->system, *relptr);
1695    
1696          /* get derivatives of y */          /* for each relation add matching jacobian elements to our matrix output */
1697          status = relman_diff3(*relptr, &vfilter, derivatives, variables, &count, enginedata->safeeval);          status = relman_diff3(*relptr, &vfilter, derivatives, variables, &count, enginedata->safeeval);
1698          if(status){          if(status){
1699              CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);              CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
# Line 1681  void integrator_ida_write_jacobian(Integ Line 1701  void integrator_ida_write_jacobian(Integ
1701              break;              break;
1702          }          }
1703    
1704          if(type==II_WRITE_YDOT){          for(j=0; j<count; ++j){
1705              for(j=0; j<count; ++j){              if(type1==II_WRITE_Y){
1706                  if(var_deriv(variables[j])){                  CONSOLE_DEBUG("looking at rel=%d, var_sindex=%d, val=%f",i,var_sindex(variables[j]),derivatives[j]);
1707                      fprintf(f, "%d %d %10.3g\n", i, integrator_ida_diffindex(blsys,variables[j]), derivatives[j]);                  fprintf(f, "%d %d %10.3g\n", i, var_sindex(variables[j]), derivatives[j]);
1708                  }              }else{
1709              }                  varname = var_make_name(sys->system,variables[j]);
1710          }else if(type==II_WRITE_Y){                  CONSOLE_DEBUG("looking at rel=%d, var_sindex=%d ('%s'), val=%f",i,integrator_ida_diffindex(sys,variables[j]),varname,derivatives[j]);
1711              for(j=0; j<count; ++j){                  fprintf(f, "%d %d %10.3g\n", i, integrator_ida_diffindex(sys,variables[j]), derivatives[j]);
1712                  if(!var_deriv(variables[j])){                  ASC_FREE(varname);
                     fprintf(f, "%d %d %10.3g\n", i, var_sindex(variables[j]), derivatives[j]);  
                 }  
1713              }              }
1714          }          }
1715      }      }
1716      ASC_FREE(variables);      ASC_FREE(variables);
1717      ASC_FREE(derivatives);      ASC_FREE(derivatives);
1718    
1719        CONSOLE_DEBUG("... DONE writing sparse Jacobian to file.d");
1720    
1721        return 0;
1722  }  }
1723    
1724  /**  /**
1725      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
1726      of debugging.      of debugging.
1727  */  */
1728  void integrator_ida_write_incidence(IntegratorSystem *blsys){  void integrator_ida_write_incidence(IntegratorSystem *sys){
1729      int i, j;      int i, j;
1730      struct rel_relation **relptr;      struct rel_relation **relptr;
1731      IntegratorIdaData *enginedata = blsys->enginedata;      IntegratorIdaData *enginedata = sys->enginedata;
1732      double *derivatives;      double *derivatives;
1733      struct var_variable **variables;      struct var_variable **variables;
1734      int count, status;      int count, status;
# Line 1717  void integrator_ida_write_incidence(Inte Line 1739  void integrator_ida_write_incidence(Inte
1739          return;          return;
1740      }      }
1741    
1742      variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);      variables = ASC_NEW_ARRAY(struct var_variable *, sys->n_y * 2);
1743      derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);      derivatives = ASC_NEW_ARRAY(double, sys->n_y * 2);
1744    
1745      CONSOLE_DEBUG("Outputting incidence information to console...");      CONSOLE_DEBUG("Outputting incidence information to console...");
1746    
# Line 1726  void integrator_ida_write_incidence(Inte Line 1748  void integrator_ida_write_incidence(Inte
1748              i< enginedata->nrels && relptr != NULL;              i< enginedata->nrels && relptr != NULL;
1749              ++i, ++relptr              ++i, ++relptr
1750      ){      ){
1751          relname = rel_make_name(blsys->system, *relptr);          relname = rel_make_name(sys->system, *relptr);
1752    
1753          /* get derivatives for this particular relation */          /* get derivatives for this particular relation */
1754          status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);          status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
# Line 1741  void integrator_ida_write_incidence(Inte Line 1763  void integrator_ida_write_incidence(Inte
1763    
1764          for(j=0; j<count; ++j){          for(j=0; j<count; ++j){
1765              if(var_deriv(variables[j])){              if(var_deriv(variables[j])){
1766                  fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(blsys,variables[j]));                  fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(sys,variables[j]));
1767              }else{              }else{
1768                  fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));                  fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
1769              }              }
# Line 1866  void integrator_ida_error(int error_code Line 1888  void integrator_ida_error(int error_code
1888          , const char *module, const char *function          , const char *module, const char *function
1889          , char *msg, void *eh_data          , char *msg, void *eh_data
1890  ){  ){
1891      IntegratorSystem *blsys;      IntegratorSystem *sys;
1892      error_severity_t sev;      error_severity_t sev;
1893    
1894      /* cast back the IntegratorSystem, just in case we need it */      /* cast back the IntegratorSystem, just in case we need it */
1895      blsys = (IntegratorSystem *)eh_data;      sys = (IntegratorSystem *)eh_data;
1896    
1897      /* severity depends on the sign of the error_code value */      /* severity depends on the sign of the error_code value */
1898      if(error_code <= 0){      if(error_code <= 0){

Legend:
Removed from v.1305  
changed lines
  Added in v.1306

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22