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

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

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

revision 995 by johnpye, Fri Dec 22 14:28:40 2006 UTC revision 996 by johnpye, Sat Dec 23 06:53:04 2006 UTC
# Line 684  int integrator_ida_solve( Line 684  int integrator_ida_solve(
684    
685          CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);          CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
686    
687          /* correct initial values, given derivatives */          /* catch SIGFPE if desired to */
688            if(enginedata->safeeval){
689                Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
690            }else{
691    #ifdef FEX_DEBUG
692                ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
693    #endif
694                Asc_SignalHandlerPushDefault(SIGFPE);
695            }
696            
697                /* correct initial values, given derivatives */
698  # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3  # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
699          /* note the new API from version 2.3 and onwards */              /* note the new API from version 2.3 and onwards */
700          flag = IDACalcIC(ida_mem, icopt, tout1);              flag = IDACalcIC(ida_mem, icopt, tout1);
701  # else  # else
702          flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);              flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
703  # endif  # endif
704    
705          switch(flag){              switch(flag){
706              case IDA_SUCCESS:                  case IDA_SUCCESS:
707                  CONSOLE_DEBUG("Initial conditions solved OK");                      CONSOLE_DEBUG("Initial conditions solved OK");
708                  break;                      break;
709    
710              case IDA_LSETUP_FAIL:                  case IDA_LSETUP_FAIL:
711              case IDA_LINIT_FAIL:                  case IDA_LINIT_FAIL:
712              case IDA_LSOLVE_FAIL:                  case IDA_LSOLVE_FAIL:
713              case IDA_NO_RECOVERY:                  case IDA_NO_RECOVERY:
714                  flag1 = -999;                      flag1 = -999;
715                  flag = (flagfn)(ida_mem,&flag1);                      flag = (flagfn)(ida_mem,&flag1);
716                  if(flag){                      if(flag){
717                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);                          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
718                            return 0;
719                        }
720                        ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
721                        return 0;
722    
723                    default:
724                        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
725                      return 0;                      return 0;
726                  }              }
727                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);  
728                  return 0;          if(enginedata->safeeval){
729                Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
730              default:          }else{
731                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");              Asc_SignalHandlerPopDefault(SIGFPE);
                 return 0;  
732          }          }
733    
734      }      }
735    
736      /* optionally, specify ROO-FINDING PROBLEM */      /* optionally, specify ROO-FINDING PROBLEM */
# Line 1007  int integrator_ida_djex(long int Neq, re Line 1024  int integrator_ida_djex(long int Neq, re
1024          /* 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) */
1025          for(j=0; j < count; ++j){          for(j=0; j < count; ++j){
1026              var_yindex = blsys->y_id[variables[j]];              var_yindex = blsys->y_id[variables[j]];
 #ifndef __WIN32__  
1027              /* the SUNDIALS headers seem not to store 'N' on Windows */              /* the SUNDIALS headers seem not to store 'N' on Windows */
1028              ASC_ASSERT_RANGE(var_yindex, -Jac->N, Jac->N);              ASC_ASSERT_RANGE(var_yindex, -blsys->n_y, blsys->n_y);
1029  #endif  
1030              if(var_yindex >= 0){              if(var_yindex >= 0){
1031                  asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);                  asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
1032                  DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];                  DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
Legend:
Removed from v.995  
changed lines
  Added in v.996
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22