chem/KPP/kpp/kpp-2.1/int/rosenbrock.c

   1
   2  #define MAX(a,b) ( ((a) >= (b)) ?(a):(b)  )
   3  #define MIN(b,c) ( ((b) <  (c)) ?(b):(c)  )
   4  #define ABS(x)   ( ((x) >=  0 ) ?(x):(-x) )
   5  #define SQRT(d)  ( pow((d),0.5)  )
   6  #define SIGN(x)  ( ((x) >=  0 ) ?[0]:(-1) )
   7
   8 /*~~> Numerical constants */
   9  #define  ZERO     (KPP_REAL)0.0
  10  #define  ONE      (KPP_REAL)1.0
  11  #define  HALF     (KPP_REAL)0.5
  12  #define  DeltaMin (KPP_REAL)1.0e-6
  13
  14 /*~~~> Collect statistics: global variables */
  15  int Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng;
  16
  17
  18 /*~~~> Function headers */
  19  void FunTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []);
  20  void JacTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []) ;
  21  int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend,
  22      KPP_REAL AbsTol[], KPP_REAL RelTol[],
  23      void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []),
  24      void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []),
  25      KPP_REAL RPAR[], int IPAR[]);
  26  int RosenbrockIntegrator(
  27      KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend ,
  28      KPP_REAL  AbsTol[], KPP_REAL  RelTol[],
  29      void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []),
  30      void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []),
  31      int ros_S,
  32      KPP_REAL ros_M[], KPP_REAL ros_E[],
  33      KPP_REAL ros_A[], KPP_REAL ros_C[],
  34      KPP_REAL ros_Alpha[],KPP_REAL  ros_Gamma[],
  35      KPP_REAL ros_ELO, char ros_NewF[],
  36      char Autonomous, char VectorTol, int Max_no_steps,
  37      KPP_REAL Roundoff, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart,
  38      KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe,
  39      KPP_REAL *Texit, KPP_REAL *Hexit );
  40  char ros_PrepareMatrix (
  41      KPP_REAL* H,
  42      int Direction,  KPP_REAL gam, KPP_REAL Jac0[],
  43      KPP_REAL Ghimj[], int Pivot[] );
  44  KPP_REAL ros_ErrorNorm (
  45      KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[],
  46      KPP_REAL AbsTol[], KPP_REAL RelTol[],
  47      char VectorTol );
  48  int  ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H);
  49  void ros_FunTimeDerivative (
  50      KPP_REAL T, KPP_REAL Roundoff,
  51      KPP_REAL Y[], KPP_REAL Fcn0[],
  52      void ode_Fun(KPP_REAL, KPP_REAL [], KPP_REAL []),
  53      KPP_REAL dFdT[] );
  54  void Fun( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] );
  55  void Jac_SP( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] );
  56  void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] );
  57  void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] );
  58  void DecompTemplate( KPP_REAL A[], int Pivot[], int* ising );
  59  void SolveTemplate( KPP_REAL A[], int Pivot[], KPP_REAL b[] );
  60  void WCOPY(int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY);
  61  void WAXPY(int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], int incY );
  62  void WSCAL(int N, KPP_REAL Alpha, KPP_REAL X[], int incX);
  63  KPP_REAL WLAMCH( char C );
  64  void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
  65              KPP_REAL ros_M[], KPP_REAL ros_E[],
  66              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
  67              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name );
  68  void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
  69              KPP_REAL ros_M[], KPP_REAL ros_E[],
  70              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
  71              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name );
  72  void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
  73              KPP_REAL ros_M[], KPP_REAL ros_E[],
  74              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
  75              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name );
  76  void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
  77              KPP_REAL ros_M[], KPP_REAL ros_E[],
  78              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
  79              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name );
  80  void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
  81              KPP_REAL ros_M[], KPP_REAL ros_E[],
  82              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
  83              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name );
  84  int  KppDecomp( KPP_REAL A[] );
  85  void KppSolve ( KPP_REAL A[], KPP_REAL b[] );
  86  void Update_SUN();
  87  void Update_RCONST();
  88
  89 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  90 void INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT )
  91 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  92 {
  93    static KPP_REAL  RPAR[20];
  94    static int  i, IERR, IPAR[20];
  95    static int Ns=0, Na=0, Nr=0, Ng=0;
  96
  97    for ( i = 0; i < 20; i++ ) {
  98      IPAR[i] = 0;
  99      RPAR[i] = ZERO;
 100    } /* for */
 101
 102
 103    IPAR[0] = 0;    /* non-autonomous */
 104    IPAR[1] = 1;    /* vector tolerances */
 105    RPAR[2] = STEPMIN; /* starting step */
 106    IPAR[3] = 5;    /* choice of the method */
 107
 108    IERR = Rosenbrock(VAR, TIN, TOUT,
 109            ATOL, RTOL,
 110            &FunTemplate, &JacTemplate,
 111            RPAR, IPAR);
 112
 113
 114    Ns=Ns+IPAR[12];
 115    Na=Na+IPAR[13];
 116    Nr=Nr+IPAR[14];
 117    Ng=Ng+IPAR[17];
 118    printf("\n Step=%d  Acc=%d  Rej=%d  Singular=%d\n",
 119          Ns,Na,Nr,Ng);
 120
 121
 122    if (IERR < 0)
 123      printf("\n Rosenbrock: Unsucessful step at T=%g: IERR=%d\n",
 124          TIN,IERR);
 125
 126    TIN = RPAR[10];      /* Exit time */
 127    STEPMIN = RPAR[11];  /* Last step */
 128
 129 } /* INTEGRATE */
 130
 131
 132 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 133 int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend,
 134         KPP_REAL AbsTol[], KPP_REAL RelTol[],
 135         void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []),
 136         void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []),
 137         KPP_REAL RPAR[], int IPAR[])
 138 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 139
 140     Solves the system y'=F(t,y) using a Rosenbrock method defined by:
 141
 142      G = 1/(H*gamma[0]) - ode_Jac(t0,Y0)
 143      T_i = t0 + Alpha(i)*H
 144      Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
 145      G * K_i = ode_Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
 146          gamma(i)*dF/dT(t0, Y0)
 147      Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
 148
 149     For details on Rosenbrock methods and their implementation consult:
 150       E. Hairer and G. Wanner
 151       "Solving ODEs II. Stiff and differential-algebraic problems".
 152       Springer series in computational mathematics, Springer-Verlag, 1996.
 153     The codes contained in the book inspired this implementation.
 154
 155     (C)  Adrian Sandu, August 2004
 156     Virginia Polytechnic Institute and State University
 157     Contact: sandu@cs.vt.edu
 158     This implementation is part of KPP - the Kinetic PreProcessor
 159 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 160
 161   *~~~>   INPUT ARGUMENTS:
 162
 163 -     Y(NVAR)    = vector of initial conditions (at T=Tstart)
 164 -    [Tstart,Tend]  = time range of integration
 165      (if Tstart>Tend the integration is performed backwards in time)
 166 -    RelTol, AbsTol = user precribed accuracy
 167 -    void ode_Fun( T, Y, Ydot ) = ODE function,
 168                        returns Ydot = Y' = F(T,Y)
 169 -    void ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
 170                        returns Jcb = dF/dY
 171 -    IPAR(1:10)    = int inputs parameters
 172 -    RPAR(1:10)    = real inputs parameters
 173 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 174
 175   *~~~>     OUTPUT ARGUMENTS:
 176
 177 -    Y(NVAR)    -> vector of final states (at T->Tend)
 178 -    IPAR(11:20)   -> int output parameters
 179 -    RPAR(11:20)   -> real output parameters
 180 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 181
 182   *~~~>    RETURN VALUE (int):
 183
 184 -    IERR       -> job status upon return
 185        - succes (positive value) or failure (negative value) -
 186            =  1 : Success
 187            = -1 : Improper value for maximal no of steps
 188            = -2 : Selected Rosenbrock method not implemented
 189            = -3 : Hmin/Hmax/Hstart must be positive
 190            = -4 : FacMin/FacMax/FacRej must be positive
 191            = -5 : Improper tolerance values
 192            = -6 : No of steps exceeds maximum bound
 193            = -7 : Step size too small
 194            = -8 : Matrix is repeatedly singular
 195 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 196
 197   *~~~>     INPUT PARAMETERS:
 198
 199     Note: For input parameters equal to zero the default values of the
 200        corresponding variables are used.
 201
 202     IPAR[0]   = 1: F = F(y)   Independent of T (AUTONOMOUS)
 203         = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
 204     IPAR[1]   = 0: AbsTol, RelTol are NVAR-dimensional vectors
 205         = 1:  AbsTol, RelTol are scalars
 206     IPAR[2]  -> maximum number of integration steps
 207         For IPAR[2]=0) the default value of 100000 is used
 208
 209     IPAR[3]  -> selection of a particular Rosenbrock method
 210         = 0 :  default method is Rodas3
 211         = 1 :  method is  Ros2
 212         = 2 :  method is  Ros3
 213         = 3 :  method is  Ros4
 214         = 4 :  method is  Rodas3
 215         = 5:   method is  Rodas4
 216
 217     RPAR[0]  -> Hmin, lower bound for the integration step size
 218           It is strongly recommended to keep Hmin = ZERO
 219     RPAR[1]  -> Hmax, upper bound for the integration step size
 220     RPAR[2]  -> Hstart, starting value for the integration step size
 221
 222     RPAR[3]  -> FacMin, lower bound on step decrease factor (default=0.2)
 223     RPAR[4]  -> FacMin,upper bound on step increase factor (default=6)
 224     RPAR[5]  -> FacRej, step decrease factor after multiple rejections
 225             (default=0.1)
 226     RPAR[6]  -> FacSafe, by which the new step is slightly smaller
 227          than the predicted value  (default=0.9)
 228 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 229
 230   *~~~>     OUTPUT PARAMETERS:
 231
 232     Note: each call to Rosenbrock adds the corrent no. of fcn calls
 233       to previous value of IPAR[10], and similar for the other params.
 234       Set IPAR(11:20) = 0 before call to avoid this accumulation.
 235
 236     IPAR[10] = No. of function calls
 237     IPAR[11] = No. of jacobian calls
 238     IPAR[12] = No. of steps
 239     IPAR[13] = No. of accepted steps
 240     IPAR[14] = No. of rejected steps (except at the beginning)
 241     IPAR[15] = No. of LU decompositions
 242     IPAR[16] = No. of forward/backward substitutions
 243     IPAR[17] = No. of singular matrix decompositions
 244
 245     RPAR[10]  -> Texit, the time corresponding to the
 246             computed Y upon return
 247     RPAR[11]  -> Hexit, last accepted step before exit
 248     For multiple restarts, use Hexit as Hstart in the following run
 249 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
 250 {
 251
 252   /*~~~>  The method parameters    */
 253    #define Smax 6
 254    int  Method, ros_S;
 255    KPP_REAL ros_M[Smax], ros_E[Smax];
 256    KPP_REAL ros_A[Smax*(Smax-1)/2], ros_C[Smax*(Smax-1)/2];
 257    KPP_REAL ros_Alpha[Smax], ros_Gamma[Smax], ros_ELO;
 258    char ros_NewF[Smax], ros_Name[12];
 259   /*~~~>  Local variables    */
 260    int Max_no_steps, IERR, i, UplimTol;
 261    char Autonomous, VectorTol;
 262    KPP_REAL Roundoff,FacMin,FacMax,FacRej,FacSafe;
 263    KPP_REAL Hmin, Hmax, Hstart, Hexit, Texit;
 264 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
 265
 266   /*~~~>  Initialize statistics */
 267    Nfun = IPAR[10];
 268    Njac = IPAR[11];
 269    Nstp = IPAR[12];
 270    Nacc = IPAR[13];
 271    Nrej = IPAR[14];
 272    Ndec = IPAR[15];
 273    Nsol = IPAR[16];
 274    Nsng = IPAR[17];
 275
 276   /*~~~>  Autonomous or time dependent ODE. Default is time dependent. */
 277    Autonomous = !(IPAR[0] == 0);
 278
 279   /*~~~>  For Scalar tolerances (IPAR[1] != 0)  the code uses AbsTol[0] and RelTol[0]
 280 !   For Vector tolerances (IPAR[1] == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR) */
 281    if (IPAR[1] == 0) {
 282       VectorTol = 1; UplimTol  = KPP_NVAR;
 283    } else {
 284       VectorTol = 0; UplimTol  = 1;
 285    } /* end if */
 286
 287   /*~~~>   The maximum number of steps admitted */
 288    if (IPAR[2] == 0)
 289       Max_no_steps = 100000;
 290    else
 291       Max_no_steps=IPAR[2];
 292    if (Max_no_steps < 0) {
 293       printf("\n User-selected max no. of steps: IPAR[2]=%d\n",IPAR[2]);
 294       return ros_ErrorMsg(-1,Tstart,ZERO);
 295    } /* end if */
 296
 297   /*~~~>  The particular Rosenbrock method chosen */
 298    if (IPAR[3] == 0)
 299        Method = 3;
 300    else
 301        Method = IPAR[3];
 302    if ( (IPAR[3] < 1) || (IPAR[3] > 5) ){
 303       printf("\n User-selected Rosenbrock method: IPAR[3]=%d\n",IPAR[3]);
 304       return ros_ErrorMsg(-2,Tstart,ZERO);
 305    } /* end if */
 306
 307   /*~~~>  Unit Roundoff (1+Roundoff>1)   */
 308    Roundoff = WLAMCH('E');
 309
 310   /*~~~>  Lower bound on the step size: (positive value) */
 311    Hmin = RPAR[0];
 312    if (RPAR[0] < ZERO) {
 313       printf("\n User-selected Hmin: RPAR[0]=%e\n", RPAR[0]);
 314       return ros_ErrorMsg(-3,Tstart,ZERO);
 315    } /* end if */
 316   /*~~~>  Upper bound on the step size: (positive value) */
 317    if (RPAR[1] == ZERO)
 318       Hmax = ABS(Tend-Tstart);
 319    else
 320       Hmax = MIN(ABS(RPAR[1]),ABS(Tend-Tstart));
 321    if (RPAR[1] < ZERO) {
 322       printf("\n User-selected Hmax: RPAR[1]=%e\n", RPAR[1]);
 323       return ros_ErrorMsg(-3,Tstart,ZERO);
 324    } /* end if */
 325   /*~~~>  Starting step size: (positive value) */
 326    if (RPAR[2] == ZERO)
 327       Hstart = MAX(Hmin,DeltaMin);
 328    else
 329       Hstart = MIN(ABS(RPAR[2]),ABS(Tend-Tstart));
 330    if (RPAR[2] < ZERO) {
 331       printf("\n User-selected Hstart: RPAR[2]=%e\n", RPAR[2]);
 332       return ros_ErrorMsg(-3,Tstart,ZERO);
 333    } /* end if */
 334   /*~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax  */
 335    if (RPAR[3] == ZERO)
 336       FacMin = (KPP_REAL)0.2;
 337    else
 338       FacMin = RPAR[3];
 339    if (RPAR[3] < ZERO) {
 340       printf("\n User-selected FacMin: RPAR[3]=%e\n", RPAR[3]);
 341       return ros_ErrorMsg(-4,Tstart,ZERO);
 342    } /* end if */
 343    if (RPAR[4] == ZERO)
 344       FacMax = (KPP_REAL)6.0;
 345    else
 346       FacMax = RPAR[4];
 347    if (RPAR[4] < ZERO) {
 348       printf("\n User-selected FacMax: RPAR[4]=%e\n", RPAR[4]);
 349       return ros_ErrorMsg(-4,Tstart,ZERO);
 350    } /* end if */
 351   /*~~~>   FacRej: Factor to decrease step after 2 succesive rejections */
 352    if (RPAR[5] == ZERO)
 353       FacRej = (KPP_REAL)0.1;
 354    else
 355       FacRej = RPAR[5];
 356    if (RPAR[5] < ZERO) {
 357       printf("\n User-selected FacRej: RPAR[5]=%e\n", RPAR[5]);
 358       return ros_ErrorMsg(-4,Tstart,ZERO);
 359    } /* end if */
 360   /*~~~>   FacSafe: Safety Factor in the computation of new step size */
 361    if (RPAR[6] == ZERO)
 362       FacSafe = (KPP_REAL)0.9;
 363    else
 364       FacSafe = RPAR[6];
 365    if (RPAR[6] < ZERO) {
 366       printf("\n User-selected FacSafe: RPAR[6]=%e\n", RPAR[6]);
 367       return ros_ErrorMsg(-4,Tstart,ZERO);
 368    } /* end if */
 369   /*~~~>  Check if tolerances are reasonable */
 370     for (i = 0; i < UplimTol; i++) {
 371       if ( (AbsTol[i] <= ZERO)  ||  (RelTol[i] <= 10.0*Roundoff)
 372           ||  (RelTol[i] >= ONE) ) {
 373         printf("\n  AbsTol[%d] = %e\n",i,AbsTol[i]);
 374         printf("\n  RelTol[%d] = %e\n",i,RelTol[i]);
 375         return ros_ErrorMsg(-5,Tstart,ZERO);
 376       } /* end if */
 377     } /* for */
 378
 379
 380   /*~~~>   Initialize the particular Rosenbrock method */
 381    switch (Method) {
 382      case 1:
 383        Ros2(&ros_S, ros_A, ros_C, ros_M, ros_E,
 384          ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name);
 385        break;
 386      case 2:
 387        Ros3(&ros_S, ros_A, ros_C, ros_M, ros_E,
 388          ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name);
 389        break;
 390      case 3:
 391        Ros4(&ros_S, ros_A, ros_C, ros_M, ros_E,
 392          ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name);
 393        break;
 394      case 4:
 395        Rodas3(&ros_S, ros_A, ros_C, ros_M, ros_E,
 396          ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name);
 397        break;
 398      case 5:
 399        Rodas4(&ros_S, ros_A, ros_C, ros_M, ros_E,
 400          ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name);
 401        break;
 402      default:
 403        printf("\n Unknown Rosenbrock method: IPAR[3]= %d", Method);
 404        return ros_ErrorMsg(-2,Tstart,ZERO);
 405    } /* end switch */
 406
 407   /*~~~>  Rosenbrock method   */
 408    IERR = RosenbrockIntegrator( Y,Tstart,Tend,
 409         AbsTol, RelTol,
 410         ode_Fun,ode_Jac ,
 411       /*  Rosenbrock method coefficients  */
 412         ros_S, ros_M, ros_E, ros_A, ros_C,
 413         ros_Alpha, ros_Gamma, ros_ELO, ros_NewF,
 414       /*  Integration parameters */
 415         Autonomous, VectorTol, Max_no_steps,
 416         Roundoff, Hmin, Hmax, Hstart,
 417         FacMin, FacMax, FacRej, FacSafe,
 418       /* Output parameters */
 419         &Texit, &Hexit );
 420
 421
 422   /*~~~>  Collect run statistics */
 423    IPAR[10] = Nfun;
 424    IPAR[11] = Njac;
 425    IPAR[12] = Nstp;
 426    IPAR[13] = Nacc;
 427    IPAR[14] = Nrej;
 428    IPAR[15] = Ndec;
 429    IPAR[16] = Nsol;
 430    IPAR[17] = Nsng;
 431   /*~~~> Last T and H */
 432    RPAR[10] = Texit;
 433    RPAR[11] = Hexit;
 434
 435    return IERR;
 436
 437 } /* Rosenbrock */
 438
 439
 440 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 441 int RosenbrockIntegrator(
 442   /*~~~> Input: the initial condition at Tstart; Output: the solution at T */
 443      KPP_REAL Y[],
 444   /*~~~> Input: integration interval */
 445      KPP_REAL Tstart, KPP_REAL Tend ,
 446   /*~~~> Input: tolerances  */
 447      KPP_REAL  AbsTol[], KPP_REAL  RelTol[],
 448   /*~~~> Input: ode function and its Jacobian */
 449      void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []),
 450      void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []) ,
 451   /*~~~> Input: The Rosenbrock method parameters */
 452      int ros_S,
 453      KPP_REAL ros_M[], KPP_REAL ros_E[],
 454      KPP_REAL ros_A[], KPP_REAL ros_C[],
 455      KPP_REAL ros_Alpha[],KPP_REAL  ros_Gamma[],
 456      KPP_REAL ros_ELO, char ros_NewF[],
 457   /*~~~> Input: integration parameters  */
 458      char Autonomous, char VectorTol,
 459      int Max_no_steps,
 460      KPP_REAL Roundoff, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart,
 461      KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe,
 462   /*~~~> Output: time at which the solution is returned (T=Tend  if success)
 463              and last accepted step  */
 464      KPP_REAL *Texit, KPP_REAL *Hexit )
 465 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 466       Template for the implementation of a generic Rosenbrock method
 467       defined by ros_S (no of stages) and coefficients ros_{A,C,M,E,Alpha,Gamma}
 468
 469       returned value: IERR, indicator of success (if positive)
 470                                       or failure (if negative)
 471 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 472 {
 473    KPP_REAL Ynew[KPP_NVAR], Fcn0[KPP_NVAR], Fcn[KPP_NVAR],
 474       dFdT[KPP_NVAR],
 475       Jac0[KPP_LU_NONZERO], Ghimj[KPP_LU_NONZERO];
 476    KPP_REAL K[KPP_NVAR*ros_S];
 477    KPP_REAL H, T, Hnew, HC, HG, Fac, Tau;
 478    KPP_REAL Err, Yerr[KPP_NVAR];
 479    int Pivot[KPP_NVAR], Direction, ioffset, j, istage;
 480    char RejectLastH, RejectMoreH;
 481
 482 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 483
 484   /*~~~>  INITIAL PREPARATIONS  */
 485    T = Tstart;
 486    *Hexit = 0.0;
 487    H = MIN(Hstart,Hmax);
 488    if (ABS(H) <= 10.0*Roundoff)
 489         H = DeltaMin;
 490
 491    if (Tend  >=  Tstart) {
 492      Direction = +1;
 493    } else {
 494      Direction = -1;
 495    } /* end if */
 496
 497    RejectLastH=0; RejectMoreH=0;
 498
 499   /*~~~> Time loop begins below  */
 500
 501    while ( ( (Direction > 0) && ((T-Tend)+Roundoff <= ZERO) )
 502        || ( (Direction < 0) && ((Tend-T)+Roundoff <= ZERO) ) ) {
 503
 504    if ( Nstp > Max_no_steps )  {                /* Too many steps */
 505         *Texit = T;
 506         return ros_ErrorMsg(-6,T,H);
 507    }
 508    if ( ((T+0.1*H) == T) || (H <= Roundoff) ) { /* Step size too small */
 509         *Texit = T;
 510         return ros_ErrorMsg(-7,T,H);
 511    }
 512
 513   /*~~~>  Limit H if necessary to avoid going beyond Tend   */
 514    *Hexit = H;
 515    H = MIN(H,ABS(Tend-T));
 516
 517   /*~~~>   Compute the function at current time  */
 518    (*ode_Fun)(T,Y,Fcn0);
 519
 520   /*~~~>  Compute the function derivative with respect to T  */
 521    if (!Autonomous)
 522       ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, ode_Fun, dFdT );
 523
 524   /*~~~>   Compute the Jacobian at current time  */
 525    (*ode_Jac)(T,Y,Jac0);
 526
 527   /*~~~>  Repeat step calculation until current step accepted  */
 528    while (1) { /* WHILE STEP NOT ACCEPTED */
 529
 530
 531    if( ros_PrepareMatrix( &H, Direction, ros_Gamma[0],
 532           Jac0, Ghimj, Pivot) ) { /* More than 5 consecutive failed decompositions */
 533        *Texit = T;
 534        return ros_ErrorMsg(-8,T,H);
 535    }
 536
 537   /*~~~>   Compute the stages  */
 538    for (istage = 1; istage <= ros_S; istage++) {
 539
 540       /* Current istage offset. Current istage vector is K[ioffset:ioffset+KPP_NVAR-1] */
 541       ioffset = KPP_NVAR*(istage-1);
 542
 543       /* For the 1st istage the function has been computed previously */
 544       if ( istage == 1 )
 545            WCOPY(KPP_NVAR,Fcn0,1,Fcn,1);
 546       else { /* istage>1 and a new function evaluation is needed at current istage */
 547         if ( ros_NewF[istage-1] ) {
 548            WCOPY(KPP_NVAR,Y,1,Ynew,1);
 549            for (j = 1; j <= istage-1; j++)
 550              WAXPY(KPP_NVAR,ros_A[(istage-1)*(istage-2)/2+j-1],
 551                    &K[KPP_NVAR*(j-1)],1,Ynew,1);
 552            Tau = T + ros_Alpha[istage-1]*Direction*H;
 553            (*ode_Fun)(Tau,Ynew,Fcn);
 554         } /*end if ros_NewF(istage)*/
 555       } /* end if istage */
 556
 557       WCOPY(KPP_NVAR,Fcn,1,&K[ioffset],1);
 558       for (j = 1; j <= istage-1; j++) {
 559          HC = ros_C[(istage-1)*(istage-2)/2+j-1]/(Direction*H);
 560          WAXPY(KPP_NVAR,HC,&K[KPP_NVAR*(j-1)],1,&K[ioffset],1);
 561       } /* for j */
 562
 563       if ((!Autonomous) && (ros_Gamma[istage-1])) {
 564         HG = Direction*H*ros_Gamma[istage-1];
 565         WAXPY(KPP_NVAR,HG,dFdT,1,&K[ioffset],1);
 566       } /* end if !Autonomous */
 567
 568       SolveTemplate(Ghimj, Pivot, &K[ioffset]);
 569
 570    } /* for istage */
 571
 572
 573   /*~~~>  Compute the new solution   */
 574    WCOPY(KPP_NVAR,Y,1,Ynew,1);
 575    for (j=1; j<=ros_S; j++)
 576        WAXPY(KPP_NVAR,ros_M[j-1],&K[KPP_NVAR*(j-1)],1,Ynew,1);
 577
 578   /*~~~>  Compute the error estimation   */
 579    WSCAL(KPP_NVAR,ZERO,Yerr,1);
 580    for (j=1; j<=ros_S; j++)
 581        WAXPY(KPP_NVAR,ros_E[j-1],&K[KPP_NVAR*(j-1)],1,Yerr,1);
 582    Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol );
 583
 584   /*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax  */
 585    Fac  = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,ONE/ros_ELO)));
 586    Hnew = H*Fac;
 587
 588   /*~~~>  Check the error magnitude and adjust step size  */
 589    Nstp++;
 590    if ( (Err <= ONE) || (H <= Hmin) ) {    /*~~~> Accept step  */
 591       Nacc++;
 592       WCOPY(KPP_NVAR,Ynew,1,Y,1);
 593       T += Direction*H;
 594       Hnew = MAX(Hmin,MIN(Hnew,Hmax));
 595       /* No step size increase after a rejected step  */
 596       if (RejectLastH)
 597          Hnew = MIN(Hnew,H);
 598       RejectLastH = 0; RejectMoreH = 0;
 599       H = Hnew;
 600          break; /* EXIT THE LOOP: WHILE STEP NOT ACCEPTED */
 601    } else {             /*~~~> Reject step  */
 602       if (Nacc >= 1)
 603          Nrej++;
 604       if (RejectMoreH)
 605          Hnew=H*FacRej;
 606       RejectMoreH = RejectLastH; RejectLastH = 1;
 607       H = Hnew;
 608    } /* end if Err <= 1 */
 609
 610    } /* while LOOP: WHILE STEP NOT ACCEPTED */
 611
 612    } /* while: time loop */
 613
 614   /*~~~> The integration was successful */
 615    *Texit = T;
 616    return 1;
 617
 618 }  /* RosenbrockIntegrator */
 619
 620
 621
 622 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 623 KPP_REAL ros_ErrorNorm (
 624   /*~~~> Input arguments */
 625      KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[],
 626      KPP_REAL AbsTol[], KPP_REAL RelTol[],
 627      char VectorTol )
 628 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 629         Computes and returns the "scaled norm" of the error vector Yerr
 630 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 631 {
 632   /*~~~> Local variables */
 633    KPP_REAL Err, Scale, Ymax;
 634    int i;
 635
 636    Err = ZERO;
 637    for (i=0; i<KPP_NVAR; i++) {
 638         Ymax = MAX(ABS(Y[i]),ABS(Ynew[i]));
 639      if (VectorTol) {
 640        Scale = AbsTol[i]+RelTol[i]*Ymax;
 641      } else {
 642        Scale = AbsTol[0]+RelTol[0]*Ymax;
 643      } /* end if */
 644      Err = Err+(Yerr[i]*Yerr[i])/(Scale*Scale);
 645    } /* for i */
 646    Err  = SQRT(Err/(KPP_REAL)KPP_NVAR);
 647
 648    return Err;
 649
 650 } /* ros_ErrorNorm */
 651
 652
 653 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 654 void ros_FunTimeDerivative (
 655     /*~~~> Input arguments: */
 656         KPP_REAL T, KPP_REAL Roundoff,
 657         KPP_REAL Y[], KPP_REAL Fcn0[],
 658         void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []),
 659     /*~~~> Output arguments: */
 660         KPP_REAL dFdT[] )
 661 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 662     The time partial derivative of the function by finite differences
 663 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 664 {
 665   /*~~~> Local variables */
 666    KPP_REAL Delta;
 667
 668    Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T));
 669    (*ode_Fun)(T+Delta,Y,dFdT);
 670    WAXPY(KPP_NVAR,(-ONE),Fcn0,1,dFdT,1);
 671    WSCAL(KPP_NVAR,(ONE/Delta),dFdT,1);
 672
 673 }  /*  ros_FunTimeDerivative */
 674
 675
 676 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 677 char ros_PrepareMatrix (
 678        /* Inout argument: (step size is decreased when LU fails) */
 679            KPP_REAL* H,
 680        /* Input arguments: */
 681            int Direction,  KPP_REAL gam, KPP_REAL Jac0[],
 682        /* Output arguments: */
 683            KPP_REAL Ghimj[], int Pivot[] )
 684 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 685   Prepares the LHS matrix for stage calculations
 686   1.  Construct Ghimj = 1/(H*ham) - Jac0
 687       "(Gamma H) Inverse Minus Jacobian"
 688   2.  Repeat LU decomposition of Ghimj until successful.
 689        -half the step size if LU decomposition fails and retry
 690        -exit after 5 consecutive fails
 691
 692   Return value:       Singular (true=1=failed_LU or false=0=successful_LU)
 693 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 694 {
 695   /*~~~> Local variables */
 696    int i, ising, Nconsecutive;
 697    KPP_REAL ghinv;
 698
 699    Nconsecutive = 0;
 700
 701    while (1) {  /* while Singular */
 702
 703   /*~~~>    Construct Ghimj = 1/(H*ham) - Jac0 */
 704      WCOPY(KPP_LU_NONZERO,Jac0,1,Ghimj,1);
 705      WSCAL(KPP_LU_NONZERO,(-ONE),Ghimj,1);
 706      ghinv = ONE/(Direction*(*H)*gam);
 707      for (i=0; i<KPP_NVAR; i++) {
 708        Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+ghinv;
 709      } /* for i */
 710   /*~~~>    Compute LU decomposition  */
 711      DecompTemplate( Ghimj, Pivot, &ising );
 712      if (ising == 0) {
 713   /*~~~>    if successful done  */
 714         return 0;  /* Singular = false */
 715      } else { /* ising .ne. 0 */
 716   /*~~~>    if unsuccessful half the step size; if 5 consecutive fails return */
 717         Nsng++; Nconsecutive++;
 718         printf("\nWarning: LU Decomposition returned ising = %d\n",ising);
 719         if (Nconsecutive <= 5) { /* Less than 5 consecutive failed LUs */
 720           *H = (*H)*HALF;
 721         } else {                  /* More than 5 consecutive failed LUs */
 722           return 1; /* Singular = true */
 723         } /* end if  Nconsecutive */
 724      } /* end if ising */
 725
 726    } /* while Singular */
 727
 728 }  /*  ros_PrepareMatrix */
 729
 730
 731 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 732 int ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H)
 733 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 734                   Handles all error messages and returns IERR = error Code
 735 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 736 {
 737    printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
 738    printf("\nForced exit from Rosenbrock due to the following error:\n");
 739
 740    switch (Code) {
 741    case -1:
 742       printf("--> Improper value for maximal no of steps"); break;
 743    case -2:
 744       printf("--> Selected Rosenbrock method not implemented"); break;
 745    case -3:
 746       printf("--> Hmin/Hmax/Hstart must be positive"); break;
 747    case -4:
 748       printf("--> FacMin/FacMax/FacRej must be positive"); break;
 749    case -5:
 750       printf("--> Improper tolerance values"); break;
 751    case -6:
 752       printf("--> No of steps exceeds maximum bound"); break;
 753    case -7:
 754       printf("--> Step size too small (T + H/10 = T) or H < Roundoff"); break;
 755    case -8:
 756       printf("--> Matrix is repeatedly singular"); break;
 757    default:
 758       printf("Unknown Error code: %d ",Code);
 759    } /* end switch */
 760
 761    printf("\n   Time = %15.7e,  H = %15.7e",T,H);
 762    printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
 763
 764    return Code;
 765
 766 }  /* ros_ErrorMsg  */
 767
 768
 769 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 770 void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
 771            KPP_REAL ros_M[], KPP_REAL ros_E[],
 772            KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
 773            char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name )
 774 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 775              AN L-STABLE METHOD, 2 stages, order 2
 776 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 777 {
 778    double g = (KPP_REAL)1.70710678118655; /* 1.0 + 1.0/SQRT(2.0) */
 779
 780   /*~~~> Name of the method */
 781     strcpy(ros_Name, "ROS-2");
 782
 783   /*~~~> Number of stages */
 784     *ros_S = 2;
 785
 786   /*~~~> The coefficient matrices A and C are strictly lower triangular.
 787     The lower triangular (subdiagonal) elements are stored in row-wise order:
 788     A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc.
 789     The general mapping formula is:
 790         A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ]   */
 791     ros_A[0] = 1.0/g;
 792
 793   /*~~~>     C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1]  */
 794     ros_C[0] = (-2.0)/g;
 795
 796   /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
 797     or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */
 798     ros_NewF[0] = 1;
 799     ros_NewF[1] = 1;
 800
 801   /*~~~> M_i = Coefficients for new step solution */
 802     ros_M[0]= (3.0)/(2.0*g);
 803     ros_M[1]= (1.0)/(2.0*g);
 804
 805   /*~~~> E_i = Coefficients for error estimator */
 806     ros_E[0] = 1.0/(2.0*g);
 807     ros_E[1] = 1.0/(2.0*g);
 808
 809   /*~~~> ros_ELO = estimator of local order - the minimum between the
 810 !    main and the embedded scheme orders plus one */
 811     *ros_ELO = (KPP_REAL)2.0;
 812
 813   /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
 814     ros_Alpha[0] = (KPP_REAL)0.0;
 815     ros_Alpha[1] = (KPP_REAL)1.0;
 816
 817   /*~~~> Gamma_i = \sum_j  gamma_{i,j}  */
 818     ros_Gamma[0] =  g;
 819     ros_Gamma[1] = -g;
 820
 821 }  /*  Ros2 */
 822
 823
 824 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 825 void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
 826            KPP_REAL ros_M[], KPP_REAL ros_E[],
 827            KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
 828            char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name )
 829 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 830        AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
 831 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 832 {
 833   /*~~~> Name of the method */
 834    strcpy(ros_Name, "ROS-3");
 835
 836   /*~~~> Number of stages */
 837    *ros_S = 3;
 838
 839   /*~~~> The coefficient matrices A and C are strictly lower triangular.
 840     The lower triangular (subdiagonal) elements are stored in row-wise order:
 841     A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc.
 842     The general mapping formula is:
 843         A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ]   */
 844    ros_A[0]= (KPP_REAL)1.0;
 845    ros_A[1]= (KPP_REAL)1.0;
 846    ros_A[2]= (KPP_REAL)0.0;
 847
 848   /*~~~>     C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1]  */
 849    ros_C[0] = (KPP_REAL)(-1.0156171083877702091975600115545);
 850    ros_C[1] = (KPP_REAL)4.0759956452537699824805835358067;
 851    ros_C[2] = (KPP_REAL)9.2076794298330791242156818474003;
 852
 853   /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
 854     or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */
 855    ros_NewF[0] = 1;
 856    ros_NewF[1] = 1;
 857    ros_NewF[2] = 0;
 858
 859   /*~~~> M_i = Coefficients for new step solution */
 860    ros_M[0] = (KPP_REAL)1.0;
 861    ros_M[1] = (KPP_REAL)6.1697947043828245592553615689730;
 862    ros_M[2] = (KPP_REAL)(-0.4277225654321857332623837380651);
 863
 864   /*~~~> E_i = Coefficients for error estimator */
 865    ros_E[0] = (KPP_REAL)0.5;
 866    ros_E[1] = (KPP_REAL)(-2.9079558716805469821718236208017);
 867    ros_E[2] = (KPP_REAL)0.2235406989781156962736090927619;
 868
 869   /*~~~> ros_ELO = estimator of local order - the minimum between the
 870 !    main and the embedded scheme orders plus 1 */
 871    *ros_ELO = (KPP_REAL)3.0;
 872
 873   /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
 874    ros_Alpha[0]= (KPP_REAL)0.0;
 875    ros_Alpha[1]= (KPP_REAL)0.43586652150845899941601945119356;
 876    ros_Alpha[2]= (KPP_REAL)0.43586652150845899941601945119356;
 877
 878   /*~~~> Gamma_i = \sum_j  gamma_{i,j}  */
 879    ros_Gamma[0]= (KPP_REAL)0.43586652150845899941601945119356;
 880    ros_Gamma[1]= (KPP_REAL)0.24291996454816804366592249683314;
 881    ros_Gamma[2]= (KPP_REAL)2.1851380027664058511513169485832;
 882
 883 }  /*  Ros3 */
 884
 885 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 886
 887
 888 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 889 void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
 890            KPP_REAL ros_M[], KPP_REAL ros_E[],
 891            KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
 892            char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name )
 893 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 894      L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
 895      L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
 896
 897       E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
 898       EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
 899       SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
 900       SPRINGER-VERLAG (1990)
 901 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 902 {
 903   /*~~~> Name of the method */
 904    strcpy(ros_Name, "ROS-4");
 905
 906   /*~~~> Number of stages */
 907    *ros_S = 4;
 908
 909   /*~~~> The coefficient matrices A and C are strictly lower triangular.
 910     The lower triangular (subdiagonal) elements are stored in row-wise order:
 911     A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc.
 912     The general mapping formula is:
 913        A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ]  */
 914    ros_A[0] = (KPP_REAL)0.2000000000000000e+01;
 915    ros_A[1] = (KPP_REAL)0.1867943637803922e+01;
 916    ros_A[2] = (KPP_REAL)0.2344449711399156;
 917    ros_A[3] = ros_A[1];
 918    ros_A[4] = ros_A[2];
 919    ros_A[5] = (KPP_REAL)0.0;
 920
 921   /*~~~>     C(i,j) = (KPP_REAL)ros_C( (i-1)*(i-2)/2 + j )  */
 922    ros_C[0] = (KPP_REAL)(-0.7137615036412310e+01);
 923    ros_C[1] = (KPP_REAL)( 0.2580708087951457e+01);
 924    ros_C[2] = (KPP_REAL)( 0.6515950076447975);
 925    ros_C[3] = (KPP_REAL)(-0.2137148994382534e+01);
 926    ros_C[4] = (KPP_REAL)(-0.3214669691237626);
 927    ros_C[5] = (KPP_REAL)(-0.6949742501781779);
 928
 929   /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
 930     or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */
 931    ros_NewF[0]  = 1;
 932    ros_NewF[1]  = 1;
 933    ros_NewF[2]  = 1;
 934    ros_NewF[3]  = 0;
 935
 936   /*~~~> M_i = Coefficients for new step solution */
 937    ros_M[0] = (KPP_REAL)0.2255570073418735e+01;
 938    ros_M[1] = (KPP_REAL)0.2870493262186792;
 939    ros_M[2] = (KPP_REAL)0.4353179431840180;
 940    ros_M[3] = (KPP_REAL)0.1093502252409163e+01;
 941
 942   /*~~~> E_i  = Coefficients for error estimator */
 943    ros_E[0] = (KPP_REAL)(-0.2815431932141155);
 944    ros_E[1] = (KPP_REAL)(-0.7276199124938920e-01);
 945    ros_E[2] = (KPP_REAL)(-0.1082196201495311);
 946    ros_E[3] = (KPP_REAL)(-0.1093502252409163e+01);
 947
 948   /*~~~> ros_ELO  = estimator of local order - the minimum between the
 949 !    main and the embedded scheme orders plus 1 */
 950    *ros_ELO  = (KPP_REAL)4.0;
 951
 952   /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
 953    ros_Alpha[0] = (KPP_REAL)0.0;
 954    ros_Alpha[1] = (KPP_REAL)0.1145640000000000e+01;
 955    ros_Alpha[2] = (KPP_REAL)0.6552168638155900;
 956    ros_Alpha[3] = (KPP_REAL)ros_Alpha[2];
 957
 958   /*~~~> Gamma_i = \sum_j  gamma_{i,j}  */
 959    ros_Gamma[0] = (KPP_REAL)( 0.5728200000000000);
 960    ros_Gamma[1] = (KPP_REAL)(-0.1769193891319233e+01);
 961    ros_Gamma[2] = (KPP_REAL)( 0.7592633437920482);
 962    ros_Gamma[3] = (KPP_REAL)(-0.1049021087100450);
 963
 964 }  /*  Ros4 */
 965
 966 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 967 void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
 968              KPP_REAL ros_M[], KPP_REAL ros_E[],
 969              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
 970              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name )
 971 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 972   --- A STIFFLY-STABLE METHOD, 4 stages, order 3
 973 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 974 {
 975   /*~~~> Name of the method */
 976    strcpy(ros_Name, "RODAS-3");
 977
 978   /*~~~> Number of stages */
 979    *ros_S = 4;
 980
 981   /*~~~> The coefficient matrices A and C are strictly lower triangular.
 982     The lower triangular (subdiagonal) elements are stored in row-wise order:
 983     A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc.
 984     The general mapping formula is:
 985         A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ]  */
 986    ros_A[0] = (KPP_REAL)0.0;
 987    ros_A[1] = (KPP_REAL)2.0;
 988    ros_A[2] = (KPP_REAL)0.0;
 989    ros_A[3] = (KPP_REAL)2.0;
 990    ros_A[4] = (KPP_REAL)0.0;
 991    ros_A[5] = (KPP_REAL)1.0;
 992
 993   /*~~~>     C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1]  */
 994    ros_C[0] = (KPP_REAL)4.0;
 995    ros_C[1] = (KPP_REAL)1.0;
 996    ros_C[2] = (KPP_REAL)(-1.0);
 997    ros_C[3] = (KPP_REAL)1.0;
 998    ros_C[4] = (KPP_REAL)(-1.0);
 999    ros_C[5] = (KPP_REAL)(-2.66666666666667); /* -8/3 */
1000
1001   /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1002     or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */
1003    ros_NewF[0]  = 1;
1004    ros_NewF[1]  = 0;
1005    ros_NewF[2]  = 1;
1006    ros_NewF[3]  = 1;
1007
1008   /*~~~> M_i = Coefficients for new step solution */
1009    ros_M[0] = (KPP_REAL)2.0;
1010    ros_M[1] = (KPP_REAL)0.0;
1011    ros_M[2] = (KPP_REAL)1.0;
1012    ros_M[3] = (KPP_REAL)1.0;
1013
1014   /*~~~> E_i  = Coefficients for error estimator */
1015    ros_E[0] = (KPP_REAL)0.0;
1016    ros_E[1] = (KPP_REAL)0.0;
1017    ros_E[2] = (KPP_REAL)0.0;
1018    ros_E[3] = (KPP_REAL)1.0;
1019
1020   /*~~~> ros_ELO  = estimator of local order - the minimum between the
1021 !    main and the embedded scheme orders plus 1 */
1022    *ros_ELO  = (KPP_REAL)3.0;
1023
1024   /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
1025    ros_Alpha[0] = (KPP_REAL)0.0;
1026    ros_Alpha[1] = (KPP_REAL)0.0;
1027    ros_Alpha[2] = (KPP_REAL)1.0;
1028    ros_Alpha[3] = (KPP_REAL)1.0;
1029
1030   /*~~~> Gamma_i = \sum_j  gamma_{i,j}  */
1031    ros_Gamma[0] = (KPP_REAL)0.5;
1032    ros_Gamma[1] = (KPP_REAL)1.5;
1033    ros_Gamma[2] = (KPP_REAL)0.0;
1034    ros_Gamma[3] = (KPP_REAL)0.0;
1035
1036 }  /*  Rodas3 */
1037
1038 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1039 void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[],
1040              KPP_REAL ros_M[], KPP_REAL ros_E[],
1041              KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[],
1042              char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name )
1043 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1044      STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1045
1046       E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1047       EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1048       SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1049       SPRINGER-VERLAG (1996)
1050 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1051 {
1052   /*~~~> Name of the method */
1053    strcpy(ros_Name, "RODAS-4");
1054
1055   /*~~~> Number of stages */
1056     *ros_S = 6;
1057
1058   /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
1059     ros_Alpha[0] = (KPP_REAL)0.000;
1060     ros_Alpha[1] = (KPP_REAL)0.386;
1061     ros_Alpha[2] = (KPP_REAL)0.210;
1062     ros_Alpha[3] = (KPP_REAL)0.630;
1063     ros_Alpha[4] = (KPP_REAL)1.000;
1064     ros_Alpha[5] = (KPP_REAL)1.000;
1065
1066   /*~~~> Gamma_i = \sum_j  gamma_{i,j}  */
1067     ros_Gamma[0] = (KPP_REAL)0.2500000000000000;
1068     ros_Gamma[1] = (KPP_REAL)(-0.1043000000000000);
1069     ros_Gamma[2] = (KPP_REAL)0.1035000000000000;
1070     ros_Gamma[3] = (KPP_REAL)(-0.3620000000000023e-01);
1071     ros_Gamma[4] = (KPP_REAL)0.0;
1072     ros_Gamma[5] = (KPP_REAL)0.0;
1073
1074   /*~~~> The coefficient matrices A and C are strictly lower triangular.
1075     The lower triangular (subdiagonal) elements are stored in row-wise order:
1076     A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc.
1077     The general mapping formula is:  A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ]  */
1078     ros_A[0]  = (KPP_REAL)0.1544000000000000e+01;
1079     ros_A[1]  = (KPP_REAL)0.9466785280815826;
1080     ros_A[2]  = (KPP_REAL)0.2557011698983284;
1081     ros_A[3]  = (KPP_REAL)0.3314825187068521e+01;
1082     ros_A[4]  = (KPP_REAL)0.2896124015972201e+01;
1083     ros_A[5]  = (KPP_REAL)0.9986419139977817;
1084     ros_A[6]  = (KPP_REAL)0.1221224509226641e+01;
1085     ros_A[7]  = (KPP_REAL)0.6019134481288629e+01;
1086     ros_A[8]  = (KPP_REAL)0.1253708332932087e+02;
1087     ros_A[9]  = (KPP_REAL)(-0.6878860361058950);
1088     ros_A[10] =  ros_A[6];
1089     ros_A[11] =  ros_A[7];
1090     ros_A[12] =  ros_A[8];
1091     ros_A[13] =  ros_A[9];
1092     ros_A[14] =  (KPP_REAL)1.0;
1093
1094   /*~~~>     C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1]  */
1095     ros_C[0]  = (KPP_REAL)(-0.5668800000000000e+01);
1096     ros_C[1]  = (KPP_REAL)(-0.2430093356833875e+01);
1097     ros_C[2]  = (KPP_REAL)(-0.2063599157091915);
1098     ros_C[3]  = (KPP_REAL)(-0.1073529058151375);
1099     ros_C[4]  = (KPP_REAL)(-0.9594562251023355e+01);
1100     ros_C[5]  = (KPP_REAL)(-0.2047028614809616e+02);
1101     ros_C[6]  = (KPP_REAL)( 0.7496443313967647e+01);
1102     ros_C[7]  = (KPP_REAL)(-0.1024680431464352e+02);
1103     ros_C[8]  = (KPP_REAL)(-0.3399990352819905e+02);
1104     ros_C[9]  = (KPP_REAL)( 0.1170890893206160e+02);
1105     ros_C[10] = (KPP_REAL)( 0.8083246795921522e+01);
1106     ros_C[11] = (KPP_REAL)(-0.7981132988064893e+01);
1107     ros_C[12] = (KPP_REAL)(-0.3152159432874371e+02);
1108     ros_C[13] = (KPP_REAL)( 0.1631930543123136e+02);
1109     ros_C[14] = (KPP_REAL)(-0.6058818238834054e+01);
1110
1111   /*~~~> M_i  = Coefficients for new step solution */
1112     ros_M[0] = ros_A[6];
1113     ros_M[1] = ros_A[7];
1114     ros_M[2] = ros_A[8];
1115     ros_M[3] = ros_A[9];
1116     ros_M[4] = (KPP_REAL)1.0;
1117     ros_M[5] = (KPP_REAL)1.0;
1118
1119   /*~~~> E_i  = Coefficients for error estimator */
1120     ros_E[0] = (KPP_REAL)0.0;
1121     ros_E[1] = (KPP_REAL)0.0;
1122     ros_E[2] = (KPP_REAL)0.0;
1123     ros_E[3] = (KPP_REAL)0.0;
1124     ros_E[4] = (KPP_REAL)0.0;
1125     ros_E[5] = (KPP_REAL)1.0;
1126
1127   /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1128     or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */
1129     ros_NewF[0] = 1;
1130     ros_NewF[1] = 1;
1131     ros_NewF[2] = 1;
1132     ros_NewF[3] = 1;
1133     ros_NewF[4] = 1;
1134     ros_NewF[5] = 1;
1135
1136   /*~~~> ros_ELO  = estimator of local order - the minimum between the
1137 !    main and the embedded scheme orders plus 1 */
1138     *ros_ELO = (KPP_REAL)4.0;
1139
1140 }  /*  Rodas4 */
1141
1142
1143
1144 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1145 void DecompTemplate( KPP_REAL A[], int Pivot[], int* ising )
1146 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1147         Template for the LU decomposition
1148 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1149 {
1150    *ising = KppDecomp ( A );
1151   /*~~~> Note: for a full matrix use Lapack:
1152       DGETRF( KPP_NVAR, KPP_NVAR, A, KPP_NVAR, Pivot, ising ) */
1153
1154    Ndec++;
1155
1156 }  /*  DecompTemplate */
1157
1158 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1159  void SolveTemplate( KPP_REAL A[], int Pivot[], KPP_REAL b[] )
1160 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1161      Template for the forward/backward substitution (using pre-computed LU decomposition)
1162 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1163 {
1164    KppSolve( A, b );
1165   /*~~~> Note: for a full matrix use Lapack:
1166       NRHS = 1
1167       DGETRS( 'N', KPP_NVAR , NRHS, A, KPP_NVAR, Pivot, b, KPP_NVAR, INFO ) */
1168
1169    Nsol++;
1170
1171 }  /*  SolveTemplate */
1172
1173
1174 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1175 void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] )
1176 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1177     Template for the ODE function call.
1178     Updates the rate coefficients (and possibly the fixed species) at each call
1179 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1180 {
1181    KPP_REAL Told;
1182
1183    Told = TIME;
1184    TIME = T;
1185    Update_SUN();
1186    Update_RCONST();
1187    Fun( Y, FIX, RCONST, Ydot );
1188    TIME = Told;
1189
1190    Nfun++;
1191
1192 }  /*  FunTemplate */
1193
1194
1195 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1196 void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] )
1197 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1198     Template for the ODE Jacobian call.
1199     Updates the rate coefficients (and possibly the fixed species) at each call
1200 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1201 {
1202   /*~~~> Local variables */
1203    KPP_REAL Told;
1204
1205    Told = TIME;
1206    TIME = T ;
1207    Update_SUN();
1208    Update_RCONST();
1209    Jac_SP( Y, FIX, RCONST, Jcb );
1210    TIME = Told;
1211
1212    Njac++;
1213
1214 } /* JacTemplate   */
1215
1216