Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/runge_kutta.c

Last change on this file was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk
File size: 75.4 KB
Line 
1 #define MAX(a,b) ( ((a) >= (b)) ?(a):(b)  )
2 #define MIN(b,c) ( ((b) <  (c)) ?(b):(c)  )
3 #define ABS(x)   ( ((x) >=  0 ) ?(x):(-x) )
4 #define SQRT(d)  ( pow((d),0.5)  )
5 /* SIGN transfer function */
6 #define SIGN(x,y) (((y) >=  0 ) ?(ABS(x)):(-ABS(x)) )
7
8/*~~> Numerical constants */
9 #define  ZERO     (KPP_REAL)0.0
10 #define  ONE      (KPP_REAL)1.0
11
12/*~~~> Statistics on the work performed by the Runge-Kutta method */
13 #define Nfun   0
14 #define Njac   1
15 #define Nstp   2
16 #define Nacc   3
17 #define Nrej   4
18 #define Ndec   5
19 #define Nsol   6
20 #define Nsng   7
21 #define Ntexit 0
22 #define Nhacc  1
23 #define Nhnew  2
24
25/*~~~> Runge-Kutta method parameters */
26 #define RKmax 3
27
28 #define R2A 1
29 #define R1A 2
30 #define L3C 3
31 #define GAU 4
32 #define L3A 5
33
34 int rkMethod,
35     SdirkError;
36 KPP_REAL rkT[RKmax][RKmax],
37        rkTinv[RKmax][RKmax],
38        rkTinvAinv[RKmax][RKmax],
39        rkAinvT[RKmax][RKmax],
40        rkA[RKmax+1][RKmax+1],
41        rkB[RKmax+1],
42        rkC[RKmax+1],
43        rkD[RKmax+1],
44        rkE[RKmax+1],
45        rkBgam[RKmax+2],
46        rkBhat[RKmax+2],
47        rkTheta[RKmax+1],
48        rkF[RKmax+2],
49        rkGamma,
50        rkAlpha,
51        rkBeta,
52        rkELO;
53/*~~~> Function headers */
54// void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT, int ICNTRL_U[], KPP_REAL RCNTRL_U[],
55//              int ISTATUS_U[], KPP_REAL RSTATUS_U[], int IERR_U);
56 void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT); 
57 void RungeKutta(int N, KPP_REAL T, KPP_REAL Tend, KPP_REAL Y[],
58               KPP_REAL RelTol[], KPP_REAL AbsTol[], KPP_REAL RCNTRL[],
59               int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[], int* IERR);
60 void RK_Integrator(int N, KPP_REAL*  T, KPP_REAL  Tend, KPP_REAL Y[],
61                   KPP_REAL AbsTol[], KPP_REAL RelTol[], int ITOL,
62                   int ISTATUS[], KPP_REAL RSTATUS[], KPP_REAL Hmin,
63                   KPP_REAL Hmax, KPP_REAL Hstart, KPP_REAL Roundoff,
64                   int Max_no_steps, int NewtonMaxit, int StartNewton,
65                   int Gustafsson, KPP_REAL ThetaMin, KPP_REAL NewtonTol,
66                   KPP_REAL FacSafe, KPP_REAL FacMax, KPP_REAL FacMin,
67                   KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int* IERR);
68 void RK_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H, int* IERR);
69 void RK_ErrorScale(int N, int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[], 
70                KPP_REAL Y[], KPP_REAL SCAL[]);
71 /*void RK_Transform(int N, KPP_REAL Tr[][RKmax], KPP_REAL Z1[], KPP_REAL Z2[],
72                KPP_REAL Z3[], KPP_REAL W1[], KPP_REAL W2[], KPP_REAL W3[]);*/
73 void RK_Interpolate(char action[], int N, KPP_REAL H, KPP_REAL Hold,
74         KPP_REAL Z1[], KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL CONT[][RKmax]);
75 void RK_PrepareRHS(int N, KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL FO[],
76        KPP_REAL Z1[], KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL R1[], KPP_REAL R2[],
77        KPP_REAL R3[]);
78 void RK_Decomp(int N, KPP_REAL H, KPP_REAL FJAC[], KPP_REAL E1[],
79                int IP1[], KPP_REAL E2R[], KPP_REAL E2I[], 
80                int IP2[], int* ISING, int ISTATUS[]);
81 void RK_Solve(int N, KPP_REAL H, KPP_REAL E1[], int IP1[], KPP_REAL E2R[],
82        KPP_REAL E2I[], int IP2[], KPP_REAL R1[], KPP_REAL R2[], KPP_REAL R3[],
83        int ISTATUS[]);
84 void RK_ErrorEstimate(int N, KPP_REAL H, KPP_REAL T, KPP_REAL Y[],
85        KPP_REAL FO[], KPP_REAL E1[], int IP1[], KPP_REAL Z1[],
86        KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL SCAL[], KPP_REAL* Err,
87        int FirstStep, int Reject, int ISTATUS[]);
88 void Radau2A_Coefficients();
89 void Lobatto3C_Coefficients ();
90 void Gauss_Coefficients();
91 void Radau1A_Coefficients();
92 void Lobatto3A_Coefficients();
93 void FUN_CHEM(KPP_REAL T, KPP_REAL V[], KPP_REAL FCT[]);
94 void JAC_CHEM(KPP_REAL T, KPP_REAL V[], KPP_REAL JF[]);
95 KPP_REAL RK_ErrorNorm(int N, KPP_REAL SCAL[], KPP_REAL DY[]);
96 void Fun(KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]);
97 void Jac_SP(KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]);
98 void WCOPY(int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY);
99 void WADD(int N, KPP_REAL Y[], KPP_REAL Z[], KPP_REAL TMP[]);
100 void WAXPY(int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], int incY );
101 KPP_REAL WLAMCH( char C );
102 void KppSolveCmplxR(KPP_REAL JVSR[], KPP_REAL JVSI[], KPP_REAL XR[], KPP_REAL XI[]);
103 int KppDecompCmplxR(KPP_REAL *JVSR, KPP_REAL *JVSI);
104 void Set2Zero(int N, KPP_REAL A[]);
105 int KppDecomp( KPP_REAL A[] );
106 void KppSolve ( KPP_REAL A[], KPP_REAL b[] );
107 void Update_SUN();
108 void Update_RCONST();
109 void Update_PHOTO();
110 
111/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
112//void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT, int ICNTRL_U[], KPP_REAL RCNTRL_U[],
113//             int ISTATUS_U[], KPP_REAL RSTATUS_U[], int IERR_U )
114void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT )
115/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
116/* RungeKutta - Fully Implicit 3-stage Runge-Kutta methods based on:
117          * Radau-2A   quadrature (order 5)                       
118          * Radau-1A   quadrature (order 5)                             
119          * Lobatto-3C quadrature (order 4)                             
120          * Gauss      quadrature (order 6)                             
121  By default the code employs the KPP sparse linear algebra routines     
122  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)       
123                                                                         
124    (C)  Adrian Sandu, August 2005                                       
125    Virginia Polytechnic Institute and State University                 
126    Contact: sandu@cs.vt.edu                                             
127    Revised by Philipp Miehe and Adrian Sandu, May 2006 
128    F90 to C translation by Tinting Jiang and Don Jacob, July 2006               
129    This implementation is part of KPP - the Kinetic PreProcessor       
130~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
131{
132   int IERR;
133   KPP_REAL RCNTRL[20],
134          RSTATUS[20],
135          T1,
136          T2;
137   int ICNTRL[20],
138       ISTATUS[20];
139   static int Ntotal = 0; /* for printing the number of steps */
140   
141   int i;
142   for ( i = 0; i < 20; i++ ) {
143     RCNTRL[i] = ZERO;
144     ICNTRL[i] = 0;
145   } /* for */
146
147   /*~~~> fine-tune the integrator: */
148   ICNTRL[1]  = 0;    /* 0=vector tolerances, 1=scalar tolerances */
149   ICNTRL[4]  = 8;    /* Max no. of Newton iterations */
150   ICNTRL[5]  = 0;    /* Starting values for Newton are interpolated(0) or zero(1) */
151   ICNTRL[9]  = 1;    /* 0 - classic or 1 - SDIRK error estimation */
152   ICNTRL[10] = 0;    /* Gustaffson(0) or classic(1) controller */ 
153   
154//   /*~~~> if optional parameters are given, and if they are >0,
155//          then use them to overwrite default settings */     
156//   if (ICNTRL_U != NULL) {
157//      for ( i = 0; i < 20; i++) {
158//         if (ICNTRL_U[i] > 0) {
159//                ICNTRL[i] = ICNTRL_U[i];
160//         } /* end if */
161//      } /* end for */
162//   } /* end if */
163//   if (RCNTRL_U != NULL) {
164//      for (i = 0; i < 20; i++) {
165//         if (RCNTRL_U[i] > 0) {
166//              RCNTRL[i] = RCNTRL_U[i];
167//         } /* end if */
168//      } /* end for */
169//   } /* end if */
170
171   T1 = TIN;
172   /*printf("T1=%f\n", T1);*/
173   T2 = TOUT;
174   /*printf("T2=%f\n", T2);*/
175   RungeKutta(NVAR, T1, T2, VAR, RTOL, ATOL, RCNTRL,ICNTRL,RSTATUS,ISTATUS, &IERR);
176   
177   Ntotal += ISTATUS[Nstp];
178   printf("NSTEPS=%d (%d)  O3=%E ", ISTATUS[Nstp], Ntotal, VAR[ind_O3]);
179
180//   /* if optional parameters are given for output
181//    * use them to store information in them */
182//   if (ISTATUS_U != NULL) {
183//      for (i = 0; i < 20; i++) {
184//         ISTATUS_U[i] = ISTATUS[i];
185//      } /* end for */
186//    } /* end if */
187//   if (RSTATUS_U != NULL) {
188//      for (i = 0; i < 20; i++) {
189//         RSTATUS_U[i] = RSTATUS[i];
190//      } /* end for */
191//   } /* end if */
192//   /*if (IERR_U != NULL) */
193//      IERR_U = IERR;
194
195   if (IERR < 0) {
196        printf("Runge-Kutta: Unsuccessful exit at T=%f(IERR=%d)", TIN, IERR);
197   } /* end if */
198   
199} /* end INTEGRATE */
200
201/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
202
203void RungeKutta(int N, KPP_REAL T, KPP_REAL Tend, KPP_REAL Y[],
204               KPP_REAL RelTol[], KPP_REAL AbsTol[], KPP_REAL RCNTRL[],
205               int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[], int* IERR)
206/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207
208    This implementation is based on the book and the code Radau5:
209
210        E. HAIRER AND G. WANNER
211        "SOLVING ORDINARY DIFFERENTIAL EQUATIONS II.
212             STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS."
213        SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
214        SPRINGER-VERLAG (1991)
215
216        UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
217        CH-1211 GENEVE 24, SWITZERLAND
218        E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
219   
220    Methods:
221         * Radau-2A     quadrature (order 5)
222         * Radau-1A     quadrature (order 5)
223         * Lobatto-3C   quadrature (order 4)
224         * Gauss        quadrature (order 6)
225
226    (C)  Adrian Sandu, August 2005
227    Virginia Polytechnic Institute and State University
228    Contact: sandu@cs.vt.edu
229    This implementation is part of KPP - the Kinetic PreProcessor
230~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231
232  *~~~>   INPUT ARGUMENTS:
233          ----------------
234    Note: For input parameters equal to zero the default values of the
235          corresponding variables are used.
236
237     N          Dimension of the system
238     T          Initial time value
239
240     Tend       Final T value (Tend-T may be positive or negative)
241
242     Y(N)       Initial values for Y
243
244     RelTol, AbsTol   Relative and absolute error tolerances.
245           for ICNTRL[1] = 0: AbsTol, RelTol are N-dimensional vectors
246                         = 1: AbsTol, RelTol are scalars
247  *~~~>   Integer input parameters:
248
249     ICNTRL[0] = not used
250     
251     ICNTRL[1] = 0: AbsTol, RelTol are NVAR-dimensional vectors
252               = 1: AbsTol, RelTol are scalars
253
254     ICNTRL[2] = RK method selection
255             = 1: Radau-2A      (the default)
256             = 2: Radau-1A
257             = 3: Lobatto-3C
258             = 4: Gauss
259             = 5: Lobatto-3A (not yet implemented)
260
261     ICNTRL[3] -> maximum number of integration steps
262        For ICNTRL[3] = 0 the default value of 10000 is used
263
264     ICNTRL[4] -> maximum number of Newton iterations
265        For ICNTRL[4] = 0 the default value of 8 is used
266
267     ICNTRL[5] -> starting values of Newton iterations:
268        ICNTRL[5] = 0: starting values are obtained from
269                     the extrapolated collocation solution
270                     (the default)
271        ICNTRL[5] = 1: starting values are zero
272
273     ICNTRL[9] -> switch for error estimation strategy
274        ICNTRL[9] = 0: one additional stage at c = 0,
275                     see Hairer (default)
276        ICNTRL[9] = 1: two additional stages at c = 0,
277                     and SDIRK at c = 1, stiffly accurate
278
279     ICNTRL[10] -> switch for step size strategy
280        ICNTRL[10] = 0: mod. predictive controller (Gustafsson, default)
281        ICNTRL[10] = 1: classical step size control
282        the choice 1 seems to produce safer results;
283        for simple problems, the choice 2 produces
284        often slightly faster runs
285
286  *~~~>   Real input parameters:
287     RCNTRL[0]  -> Hmin, lower bound for the integration step size
288                (highly recommended to keep Hmin = ZERO, the default)
289
290     RCNTRL[1]  -> Hmax, upper bound for the integration step size
291
292     RCNTRL[2]  -> Hstart, the starting step size
293
294     RCNTRL[3]  -> FacMin, lower bound on step decrease factor
295                (default=0.2)
296
297     RCNTRL[4]  -> FacMax, upper bound on step increase factor
298                (default=6)
299
300     RCNTRL[5]  -> FacRej, step decrease factor after multiple rejections
301                (default=0.1)
302
303     RCNTRL[6]  -> FacSafe, by which the new step is slightly smaller
304                than the predicted value (default=0.9)
305
306     RCNTRL[7]  -> ThetaMin. If Newton convergence rate smaller
307                than ThetaMin the Jacobian is not recomputed;
308                (default=0.001)
309
310     RCNTRL[8]  -> NewtonTol, stopping criterion for Newton's method
311                (default=0.03)
312
313     RCNTRL[9]  -> Qmin
314
315     RCNTRL[10] -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
316                 step size is kept constant and the LU factorization
317                 reused (default Qmin=1, Qmax=1.2)     
318~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
319
320  *~~~>     OUTPUT ARGUMENTS:
321            -----------------
322     T        -> T value for which the solution has been computed
323                 (after successful return T=Tend).
324
325     Y(N)  -> Numerical solution at T
326
327  Note: each call to RungeKutta adds the current no. of fcn calls
328        to previous value of ISTATUS[0], and similar to other params.
329        Set ISTATUS[1:9] = 0 before call to avoid this accumulation.
330
331     ISTATUS[0] -> No. of function calls
332     ISTATUS[1] -> No. of Jacobian calls
333     ISTATUS[2] -> No. of steps
334     ISTATUS[3] -> No. of accepted steps
335     ISTATUS[4] -> No. of rejected steps (except at very beginning)
336     ISTATUS[5] -> No. of LU decompositions
337     ISTATUS[6] -> No. of forward/backward substitutions
338     ISTATUS[7] -> No. of singular matrix decompositions
339
340     RSTATUS[0] -> Texit, the time corresponding to the
341                 computed Y upon return
342     RSTATUS[1] -> Hexit, last accepted step before exit
343     RSTATUS[2] -> Hnew, last predicted step (not yet taken)
344                 For multiple restarts, use Hnew as Hstart
345                 in the subsequent run
346
347~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
348
349  *~~~>    RETURN VALUE (int):
350
351     IERR     -> Reports on successfulness upon return:
352        = 1 for success
353        < 0 for error (value equals error code)
354        = -1  : Improper value for maximal no of steps
355        = -2  : Improper value for maximal no of Newton iterations
356        = -3  : Hmin/Hmax/Hstart must be positive
357        = -4  : Improper values for FacMin/FacMax/FacSafe/FacRej
358        = -5  : Improper value for ThetaMin
359        = -6  : Newton stopping tolerance too small
360        = -7  : Improper values for Qmin, Qmax
361        = -8  : Tolerances are too small
362        = -9  : No of steps exceeds maximum bound
363        = -10 : Step size too small
364        = -11 : Matrix is repeatedly singular
365        = -12 : Non-convergence of Newton iterations
366        = -13 : Requested RK method not implemented
367~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
368{
369  /*~~~>  Control arguments */
370  /*printf("Starting RungeKutta\n");*/
371   int Max_no_steps;
372   KPP_REAL Hmin,
373          Hmax,
374          Hstart,
375          Qmin,
376          Qmax,
377          Roundoff,
378          ThetaMin,
379          NewtonTol,
380          FacSafe,
381          FacMin,
382          FacMax,
383          FacRej;
384  /*~~~> Local variables */
385   int NewtonMaxit,
386       ITOL,
387       i,
388       StartNewton,
389       Gustafsson; 
390
391   *IERR = 0;
392   for (i = 0; i < 20; i++) {
393        ISTATUS[i] = 0;
394        RSTATUS[i] = ZERO;
395   } /* end for */
396   
397/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
398
399  /*~~~> ICNTRL[0] - autonomous system - not used */
400  /*~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol */
401   if (ICNTRL[1] == 0) {
402        /*printf("Entering if ICNTRL[1] == 0\n");*/
403        ITOL = 1;
404   }
405   else {
406        ITOL = 0;
407   } /*end if */
408  /*~~~> Error control selection */
409   if (ICNTRL[9] == 0) {
410        SdirkError = 0;
411   }
412   else {
413        /*printf("Entering if ICNTRL[9] == 1\n");*/   
414        SdirkError = 1;
415   } /* end if */
416  /*~~~> Method Selection */
417  /*printf("Starting case ICNTRL[2]\n");*/
418   switch (ICNTRL[2]) {
419        case 0:
420        case 1: Radau2A_Coefficients(); /*printf("case 0 or 1\n");*/
421                break;
422        case 2: Lobatto3C_Coefficients(); printf("case 2\n");
423                break;
424        case 3: Gauss_Coefficients(); printf("case 3\n");
425                break;
426        case 4: Radau1A_Coefficients(); printf("case 4\n");
427                break;
428        case 5: Lobatto3A_Coefficients(); printf("case 5\n");
429                break;
430        default:
431                printf("\n ICNTRL[2]=%d\n", ICNTRL[2]);
432                RK_ErrorMsg(-13, T, ZERO, IERR);
433   } /* end switch */
434 
435 /*~~~> Max_no_steps: the maximal number of time steps */
436   if (ICNTRL[3] == 0) {
437        Max_no_steps = 200000;
438        /*printf("Max_no_steps = %d\n", Max_no_steps);*/
439   }
440   else {
441        Max_no_steps = ICNTRL[3];
442        if (Max_no_steps <= 0)
443        {
444                printf("\n ICNTRL[3]=%d\n", ICNTRL[3]);
445                RK_ErrorMsg(-1, T, ZERO, IERR);
446        } /* end if */
447   } /* end if */
448
449 /*~~~> NewtonMaxit: maximal number of Newton iterations */
450   if (ICNTRL[4] == 0)
451        NewtonMaxit = 8;
452   else 
453   {
454        NewtonMaxit = ICNTRL[4];
455        if (NewtonMaxit <= 0) {
456                printf("\n ICNTRL[4]=%d\n", ICNTRL[4]);
457                RK_ErrorMsg(-2, T, ZERO, IERR);
458        }
459        /*printf("NewtonMaxit = %d\n", NewtonMaxit);*/
460   } /* end if */
461
462 /*~~~> StartNewton: Use extrapolation for starting values of Newton iterations */
463   if (ICNTRL[5] == 0) {
464        StartNewton = 1;
465        /*printf("StartNewton = %d\n", StartNewton);*/
466   }
467   else {
468        StartNewton = 0;
469   } /* end if */
470
471 /*~~~> Gustafsson: step size controller */
472   if (ICNTRL[10] == 0) {
473        Gustafsson = 1;
474        /*printf("Gustafsson = %d\n", Gustafsson);*/
475   }
476   else {
477        Gustafsson = 0;
478   } /* end if */
479 
480 /*~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 */
481   Roundoff = WLAMCH('E');
482   /*printf("Roundoff=%f\n", Roundoff);*/
483 /*~~~> Hmin = minimal step size */
484   if (RCNTRL[0] == ZERO) {
485        /*printf("Entering if RCNTRL[0]=%f == ZERO\n", RCNTRL[0]);*/   
486        Hmin = ZERO;
487   }
488   else {
489        Hmin = MIN(ABS(RCNTRL[0]), ABS(Tend-T));
490   } /* end if */
491 /*~~~> Hmax = maximal step size */
492   if (RCNTRL[1] == ZERO) {
493        /*printf("Entering if RCNTRL[1]=%f == ZERO\n", RCNTRL[1]);*/   
494        Hmax = ABS(Tend-T);
495   }
496   else {
497        Hmax = MIN(ABS(RCNTRL[1]),ABS(Tend-T));
498   } /* end if */
499 /*~~~> Hstart = starting step size */
500   if (RCNTRL[2] == ZERO) {
501        /*printf("Entering if RCNTRL[2]=%f == ZERO\n", RCNTRL[2]);*/   
502        Hstart = ZERO;
503   }
504   else {
505        Hstart = MIN(ABS(RCNTRL[2]),ABS(Tend-T));
506   } /* end if */
507 /*~~~> FacMin: lower bound on step decrease factor */
508   if (RCNTRL[3] == ZERO) {
509        /*printf("Entering if RCNTRL[3]=%f == ZERO\n", RCNTRL[3]);*/   
510        FacMin = (KPP_REAL)0.2;
511   }
512   else {
513        FacMin = RCNTRL[3];
514   } /* end if */
515 /*~~~> FacMax: upper bound on step increase factor */
516   if (RCNTRL[4] == ZERO) {
517        /*printf("Entering if RCNTRL[4]=%f == ZERO\n", RCNTRL[4]);*/   
518        FacMax = (KPP_REAL)8.0;
519   }
520   else {
521        FacMax = RCNTRL[4];
522   } /* end if */
523 /*~~~> FacRej: step decrease factor after 2 consecutive rejections */
524   if (RCNTRL[5] == ZERO) {
525        /*printf("Entering if RCNTRL[5]=%f == ZERO\n", RCNTRL[5]);*/
526        FacRej = (KPP_REAL)0.1;
527   }
528   else {
529        FacRej = RCNTRL[5];
530   } /* end if */
531 /*~~~> FacSafe: by which the new step is slightly smaller
532                 than the predicted value */
533   if (RCNTRL[6] == ZERO) {
534        /*printf("Entering if RCNTRL[6]=%f == ZERO\n", RCNTRL[6]);*/   
535        FacSafe = (KPP_REAL)0.9;
536   }
537   else {
538        FacSafe = RCNTRL[6];
539   } /* end if */
540   if ((FacMax < ONE) || (FacMin > ONE) || (FacSafe <= 1.0e-03) || (FacSafe >= ONE)) {
541        printf("\n RCNTRL[3]=%f, RCNTRL[4]=%f, RCNTRL[5]=%f, RCNTRL[6]=%f\n",
542                RCNTRL[3], RCNTRL[4], RCNTRL[5], RCNTRL[6]);
543        RK_ErrorMsg(-4, T, ZERO, IERR);
544   } /* end if */
545
546 /*~~~> ThetaMin: decides whether the Jacobian should be recomputed */
547   if (RCNTRL[7] == ZERO)
548        ThetaMin = (KPP_REAL)1.0e-03;
549   else {
550        ThetaMin = RCNTRL[7];
551        if (ThetaMin <= (KPP_REAL)0.0 || ThetaMin >= (KPP_REAL)1.0) {
552                printf("\n RCNTRL[7]=%f\n", RCNTRL[7]);
553                RK_ErrorMsg(-5, T, ZERO, IERR);
554        }
555   } /* end if */
556 /*~~~> NewtonTol: stopping criterion for Newton's method */
557   if (RCNTRL[8] == ZERO)
558        NewtonTol = (KPP_REAL)3.0e-02;
559   else {
560        NewtonTol = RCNTRL[8];
561        if (NewtonTol <= Roundoff) {
562                printf("\n RCNTRL[8]=%f\n", RCNTRL[8]);
563                RK_ErrorMsg(-6, T, ZERO, IERR);
564        }
565   } /* end if */
566 /*~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax then step size = const. */
567   if (RCNTRL[9] == ZERO) {
568        Qmin = ONE;
569   }
570   else {
571        Qmin = RCNTRL[9];
572   } /* end if */
573   if (RCNTRL[10] == ZERO) {
574        Qmax = (KPP_REAL)1.2;
575   }
576   else {
577        Qmax = RCNTRL[10];
578   } /* end if */
579   if (Qmin > ONE || Qmax < ONE) {
580        printf("\n RCNTRL[9]=%f\n", Qmin);
581        printf("\n RCNTRL[10]=%f\n", Qmax);
582        RK_ErrorMsg(-7, T, ZERO, IERR);
583   } /* end if */
584 /*~~~> Check if tolerances are reasonable */
585   if (ITOL == 0) {
586        if ( AbsTol[0] <= ZERO || RelTol[0] <= ( ((KPP_REAL)10.0)*Roundoff) ) {
587           printf("\n AbsTol=%f\n", AbsTol[0]);
588           printf("\n RelTol=%f\n", RelTol[0]);
589           RK_ErrorMsg(-8, T, ZERO, IERR);
590        }
591   }
592   else {
593        for (i = 0; i < N; i++) {
594           if ( (AbsTol[i] <= ZERO) || (RelTol[i] <= ((KPP_REAL)10.0)*Roundoff) ) {
595                printf("\n AbsTol[%d] = %f\n", i, AbsTol[i]);
596                printf("\n AbsTol[%d] = %f\n", i, RelTol[i]);
597                RK_ErrorMsg(-8, T, ZERO, IERR);
598           } /* end if */
599        } /* end for */
600   } /* end if */
601
602 /*~~~> Parameters are wrong */
603   if (*IERR < 0)
604        return;
605
606 /*~~~> Call the core method */ 
607   RK_Integrator(N, &T, Tend, Y, AbsTol, RelTol, ITOL, ISTATUS, RSTATUS,
608                   Hmin, Hmax, Hstart, Roundoff, Max_no_steps, NewtonMaxit,
609                   StartNewton, Gustafsson, ThetaMin, NewtonTol,
610                   FacSafe, FacMax, FacMin, FacRej, Qmin, Qmax, IERR);
611  /*printf("*IERR = %d\n", *IERR);*/
612  /*printf("Ending RungeKutta\n");*/
613} /* RungeKutta */
614   
615/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
616void RK_Integrator( int N, 
617        /*~~~> Input: integration interval */       
618        KPP_REAL*  T, KPP_REAL  Tend, KPP_REAL Y[],
619        KPP_REAL AbsTol[], KPP_REAL RelTol[], int ITOL,
620  /*~~~> Input: the initial condition at T; output: the solution at Tend */     
621        int ISTATUS[], KPP_REAL RSTATUS[], KPP_REAL Hmin, KPP_REAL Hmax, 
622        KPP_REAL Hstart, KPP_REAL Roundoff, int Max_no_steps, int NewtonMaxit,
623        int StartNewton, int Gustafsson, KPP_REAL ThetaMin,
624        KPP_REAL NewtonTol, KPP_REAL FacSafe, KPP_REAL FacMax, KPP_REAL FacMin,
625        KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int* IERR)
626{ 
627   /*printf("Starting RK_Integrator\n");*/
628   KPP_REAL FJAC[LU_NONZERO],
629          E1[LU_NONZERO],
630          E2R[LU_NONZERO],
631          E2I[LU_NONZERO];
632   KPP_REAL Z1[NVAR],
633          Z2[NVAR],
634          Z3[NVAR],
635          Z4[NVAR],
636          SCAL[NVAR],
637          DZ1[NVAR],
638          DZ2[NVAR],
639          DZ3[NVAR],
640          DZ4[NVAR],
641          G[NVAR],
642          TMP[NVAR],
643          FO[NVAR];   
644   KPP_REAL CONT[NVAR][RKmax],
645          Tdirection,
646          H,
647          Hacc,
648          Hnew,
649          Hold,
650          Fac,
651          FacGus, 
652          Theta,
653          Err,
654          ErrOld,
655          NewtonRate,
656          NewtonIncrement,
657          Hratio,
658          Qnewton, 
659          NewtonPredictedErr,
660          NewtonIncrementOld,
661          ThetaSD;
662   int IP1[NVAR],
663       IP2[NVAR],
664       NewtonIter,
665       Nconsecutive,
666       i;
667   int ISING;
668   int Reject,
669       FirstStep,
670       SkipJac,
671       NewtonDone,
672       SkipLU;
673
674/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
675   
676  /*~~~>  INITIAL setting */
677   /*printf("ONE=%f Tend-*T=%f\n", ONE, Tend-*T);*/
678   Tdirection = SIGN(ONE, Tend-*T);
679   /*printf("Tdirection=%f\n", Tdirection);*/
680   /*printf("Hmin=%g Hstart=%g \n", Hmin, Hstart);*/
681   /*printf("ABS(Hmin)=%g ABS(Hstart)=%g \n", ABS(Hmin), ABS(Hstart));*/
682   H = MIN( MAX(ABS(Hmin), ABS(Hstart)), Hmax ); 
683   if (ABS(H) <= ((KPP_REAL)10.0*Roundoff)) {
684        /*printf("Entering if ABS(H)=%g <= 10.0*Roundoff=%g\n", ABS(H),(KPP_REAL)10.0*Roundoff);*/
685        H = (KPP_REAL)(1.0e-06);
686   } /* end if */
687   H = SIGN(H, Tdirection);
688   Hold = H;
689   /*printf("Hold=%f\n", Hold);*/ 
690   Reject = 0;
691   FirstStep = 1;
692   SkipJac = 0;
693   SkipLU = 0;
694   if ((*T+H*((KPP_REAL)1.0001)-Tend)*Tdirection >= ZERO) {   
695        H = Tend - *T;
696   } /* end if */
697   Nconsecutive = 0;
698   RK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL);
699   /*for(i=0; i<NVAR; i++) {
700        printf("AbsTol=%g RelTol=%g Y=%g SCAL=%g \n", AbsTol[i], RelTol[i], Y[i], SCAL[i] );
701   }*/
702   
703 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
704  /*~~~> Time loop begins */ 
705  /* while Tloop */
706Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
707      /*printf("Starting Tloop: (Tend-*T)*Tdirection - Roundoff=%g > ZERO\n", (Tend-*T)*Tdirection-Roundoff);*/
708      /*if ( Reject == 0 ) { */
709        FUN_CHEM(*T,Y,FO);
710        ISTATUS[Nfun]++;
711      /* } * end if */ 
712      if ( SkipLU == 0 ) { /* This time around skip the Jac update and LU */
713        /*~~~> Compute the Jacobian matrix */
714        /*printf("Entering if SkipLU == 0\n");*/
715        if ( SkipJac == 0 ) {
716           /*printf("Entering if SkipJac == 0\n");*/
717           JAC_CHEM(*T,Y,FJAC);
718           ISTATUS[Njac]++;
719        } /* end if */
720        /*~~~> Compute the matrices E1 and E2 and their decompositions */
721        RK_Decomp(N,H,FJAC,E1,IP1,E2R,E2I,IP2,&ISING,ISTATUS);
722        /*printf("ISING=%d\n", ISING);*/
723        if ( ISING != 0 ) {
724           /*printf("Entering if ISING != 0\n");*/
725           ISTATUS[Nsng]++;
726           Nconsecutive++;
727           if (Nconsecutive >= 5) {
728                RK_ErrorMsg(-12,*T,H,IERR);
729           }
730           H = H * ((KPP_REAL)0.5);
731           Reject = 1;
732           SkipJac = 1;
733           SkipLU = 0;
734           goto Tloop;
735        }
736        else {
737           /*printf("Entering if ISING == 0\n");*/
738           Nconsecutive = 0;
739        } /* end if */
740      } /* end if !SkipLU */
741
742      /*printf("NSTEPS=%d\n", ISTATUS[Nstp]);*/
743      ISTATUS[Nstp]++;
744      /*printf("NSTEPS=%d\n", ISTATUS[Nstp]);*/
745      if (ISTATUS[Nstp] > Max_no_steps) {
746        printf("\n Max number of time steps is = %d", Max_no_steps);
747        RK_ErrorMsg(-9,*T,H,IERR);
748      } /* end if */
749      if (((KPP_REAL)0.1)*ABS(H) <= ABS(*T)*Roundoff)
750        RK_ErrorMsg(-10,*T,H,IERR);
751 
752 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
753 /*~~~>  Loop for the simplified Newton iterations */
754 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
755
756      /*~~~> Starting values for Newton iteration */ 
757      if ( (FirstStep == 1) || (StartNewton == 0) ) {
758        /*printf("Entering if FirstStep or Not StartNewton\n");*/
759        Set2Zero(N,Z1);
760        Set2Zero(N,Z2);
761        Set2Zero(N,Z3);
762      }
763      else
764      {
765        /* Evaluate quadratic polynomial */
766        RK_Interpolate("eval",N,H,Hold,Z1,Z2,Z3,CONT);
767      } /* end if */
768
769      /*~~~> Initializations for Newton iteration */
770      NewtonDone = 0;
771      Fac = (KPP_REAL)0.5; /* Step reduction if too many iterations */
772 
773      /* for NewtonLoop */
774      for (NewtonIter = 1; NewtonIter <= NewtonMaxit; NewtonIter++)
775      {
776        /*printf("Starting NewtonIter loop: NewtonIter=%d\n", NewtonIter);*/
777        /* Prepare the right-hand side */
778        RK_PrepareRHS(N,*T,H,Y,FO,Z1,Z2,Z3,DZ1,DZ2,DZ3);
779
780        /* Solve the linear systems */
781        RK_Solve(N,H,E1,IP1,E2R,E2I,IP2,DZ1,DZ2,DZ3,ISTATUS);
782
783        NewtonIncrement = 
784                SQRT((RK_ErrorNorm(N,SCAL,DZ1) * RK_ErrorNorm(N,SCAL,DZ1)
785                    + RK_ErrorNorm(N,SCAL,DZ2) * RK_ErrorNorm(N,SCAL,DZ2) 
786                    + RK_ErrorNorm(N,SCAL,DZ3) * RK_ErrorNorm(N,SCAL,DZ3))
787                    / (KPP_REAL)3.0 );
788        /*printf( "NewtonIncrement=%g \n", NewtonIncrement );*/
789        if (NewtonIter == 1) {
790           /*printf("Entering NewtonIter=%d == 1\n", NewtonIter);*/
791           Theta = ABS(ThetaMin);
792           NewtonRate = (KPP_REAL)2.0;
793        }
794        else {
795           /*printf("Entering else from NewtonIter=%d == 1\n", NewtonIter);*/
796           Theta = NewtonIncrement / NewtonIncrementOld;
797           /*printf("Theta=%f\n", Theta);*/
798           if (Theta < (KPP_REAL)0.99) {
799                /*printf( "Entering if Theta=%g < 0.99\n", Theta );*/
800                NewtonRate = Theta / (ONE-Theta);
801           }
802           else { /* Non-convergence of Newton: Theta too large */
803              break; /* EXIT NewtonLoop */
804           } /* end if */
805           if (NewtonIter < NewtonMaxit) {
806              /*printf("Entering if NewtonIter=%d < NewtonMaxit=%d\n", NewtonIter, NewtonMaxit);*/
807              /* Predict error at the end of Newton process */
808              NewtonPredictedErr = NewtonIncrement
809                             *pow(Theta,(NewtonMaxit-NewtonIter))/(ONE-Theta);
810              /*printf("NewtonRate=%g NewtonPredictedErr=%g \n", NewtonRate, NewtonPredictedErr);*/
811              if (NewtonPredictedErr >= NewtonTol) {
812                /*printf( "Entering if NewtonPredictedErr=%g >= NewtonTol%g \n", NewtonPredictedErr, NewtonTol );*/
813                /* Non-convergence of Newton: predicted error too large */
814                Qnewton = MIN((KPP_REAL)10.0, NewtonPredictedErr/NewtonTol);
815                Fac = (KPP_REAL)0.8*pow(Qnewton,(-ONE/(1+NewtonMaxit-NewtonIter)));
816                break; /* EXIT NewtonLoop */
817              } /* end if */
818           } /* end if */
819        } /* end if */
820           
821        NewtonIncrementOld = MAX(NewtonIncrement, Roundoff);
822        /*printf("NewtonIncrementOld=%f\n", NewtonIncrementOld);*/
823        /*~~~> Update solution */
824        WAXPY(N,-ONE,DZ1,1,Z1,1);       /* Z1 <- Z1 - DZ1 */
825        WAXPY(N,-ONE,DZ2,1,Z2,1);       /* Z2 <- Z2 - DZ2 */
826        WAXPY(N,-ONE,DZ3,1,Z3,1);       /* Z3 <- Z3 - DZ3 */
827        /*for(i=0; i<N; i++)
828           printf("DZ1[%d]=%g Z1[%d]=%g\n", i, DZ1[i], i, Z1[i]);*/
829        /*~~~> Check error in Newton iterations */
830        /*printf( "NewtonRate=%g NewtonIncrement=%g NewtonTol=%g \n", NewtonRate, NewtonIncrement, NewtonTol );*/
831        NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol);
832        /*printf( "NewtonDone=%d \n", NewtonDone );*/
833        if (NewtonDone == 1)
834           break; /* Exit NewtonLoop */
835        if (NewtonIter == NewtonMaxit)
836        {
837           printf("\n Slow or no convergence in Newton Iteration:");
838           printf(" Max no. of Newton iterations reached");
839        } /* end if */
840      }/* end for NewtonLoop */
841
842      if ( NewtonDone == 0)
843      {
844        /*printf( "Entering if NewtonDone == 0\n" );*/
845        /*RK_ErrorMsg(-12,*T,H,IERR); */
846        H = Fac*H;
847        Reject  = 1;
848        SkipJac = 1;
849        SkipLU  = 0;
850        goto Tloop;
851      } /* end if */
852
853 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
854 /*~~~> SDIRK Stage */
855 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
856      if (SdirkError == 1)
857      {
858        /*printf("Entering if SdirkError=%d==1\n", SdirkError);*/
859        /*~~~> Starting values for Newton iterations */
860        for (i = 0; i < N; i++)
861        {
862           Z4[i] = Z3[i];
863        }
864        /*~~~> Prepare the loop-independent part of the right-hand side */
865        /* G = H*rkBgam(0)*FO + rkTheta(1)*Z1
866                + rkTheta(2)*Z2 + rkTheta(3)*Z3; */
867        Set2Zero(N,G);
868        if (rkMethod != L3A)
869           WAXPY(N,rkBgam[0]*H, FO,1,G,1);
870        WAXPY(N,rkTheta[0],Z1,1,G,1); 
871        WAXPY(N,rkTheta[1],Z2,1,G,1); 
872        WAXPY(N,rkTheta[2],Z3,1,G,1); 
873
874        /*~~~> Initializations for Newton iteration */
875        NewtonDone = 0;
876        Fac = (KPP_REAL)0.5; /* Step reduction factor if too many iterations */
877
878        /* for SDNewtonLoop */
879        for (NewtonIter = 1; NewtonIter <= NewtonMaxit; NewtonIter++)
880        {
881           /*printf("Starting NewtonIter loop: NewtonIter=%d\n", NewtonIter);*/
882           /*~~~> Prepare the loop-dependent part of right-hand side */
883           WADD(N,Y,Z4,TMP);    /* TMP <- Y + Z4 */
884           FUN_CHEM(*T+H,TMP,DZ4);      /* DZ4 <- Fun(Y+Z4) */
885           ISTATUS[Nfun]++;
886           /* DZ4[1,N] = (D[1, N]-Z4[1,N])*(rkGamma/H) + DZ4[1,N]; */
887           WAXPY(N,-ONE*rkGamma/H,Z4,1,DZ4,1);
888           WAXPY(N,rkGamma/H,G,1,DZ4,1);
889
890           /*~~~> Solve the linear system */
891           KppSolve(E1, DZ4);
892           /*~~~> Note: for a full matrix use Lapack:
893                DGETRS('N', 5, 1, E1, N, IP1, DZ4, 5, ISING) */
894
895           /*~~~> Check convergence of Newton iterations */
896           NewtonIncrement = RK_ErrorNorm(N,SCAL,DZ4);
897           /*printf( "NewtonIncrement=%g \n", NewtonIncrement );*/
898           if (NewtonIter == 1) {
899              /*printf("Entering NewtonIter=%d == 1\n", NewtonIter);*/
900              ThetaSD = ABS(ThetaMin);
901              NewtonRate = (KPP_REAL)2.0;
902           }
903           else {
904              /*printf("Entering else from NewtonIter=%d == 1\n", NewtonIter);*/
905              ThetaSD = NewtonIncrement / NewtonIncrementOld;
906              if (ThetaSD < (KPP_REAL)0.99) {
907                /*printf( "Entering if Theta=%g < 0.99\n", Theta );*/
908                NewtonRate = ThetaSD / (ONE-ThetaSD);
909                /* Predict error at the end of Newton process */
910                NewtonPredictedErr = NewtonIncrement
911                        *pow(ThetaSD,(NewtonMaxit-NewtonIter))/(ONE-ThetaSD);
912                /*printf("NewtonRate=%g NewtonPredictedErr=%g \n", NewtonRate, NewtonPredictedErr);*/
913                if (NewtonPredictedErr >= NewtonTol) {
914                   /*printf( "Entering if NewtonPredictedErr=%g >= NewtonTol%g \n", NewtonPredictedErr, NewtonTol );*/
915                   /* Non-convergence of Newton: predicted error too large */
916                   /* printf("\n Error too large", NewtonPredictedErr); */
917                   Qnewton = MIN((KPP_REAL)10.0,NewtonPredictedErr/NewtonTol);
918                   Fac = (KPP_REAL)0.8*pow(Qnewton,(-ONE/(1+NewtonMaxit-NewtonIter)));
919                   break; /* EXIT SDNewtonLoop */
920                } /* end if */
921              }
922              /* Non-convergence of Newton: Theta too large */
923              else {
924                /* prinf("\n Theta too large", ThetaSD); */
925                break; /* EXIT SDNewtonLoop */
926              } /* end if */
927           } /* end if */
928           NewtonIncrementOld = NewtonIncrement;
929           /*printf("NewtonIncrementOld=%f\n", NewtonIncrementOld);*/
930           /* Update solution: Z4 <-- Z4 + DZ4; */
931           WAXPY(N,ONE,DZ4,1,Z4,1);
932
933           /* Check error in Newton iterations */
934           NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol);
935           /*printf( "NewtonDone=%d \n", NewtonDone );*/
936           if (NewtonDone == 1)
937                break;  /* EXIT SDNewtonLoop */
938        } /* end for SDNewtonLoop */
939     
940        if ( NewtonDone == 0 ) {
941           /*printf( "Entering if NewtonDone == 0\n" );*/
942           H       = Fac*H;
943           Reject  = 1;
944           SkipJac = 1;
945           SkipLU  = 0;
946           goto Tloop;
947        } /* end if */
948      } /* end if */
949      /*~~~> End of implified SDIRK Newton iterations */
950
951 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
952 /*~~~> Error estimation */
953 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
954      if (SdirkError == 1) {
955        Set2Zero(N, DZ4);
956        if (rkMethod == L3A) {
957           for (i = 0; i < N; i++)
958                DZ4[i] = H*rkF[0]*FO[i];
959           if (rkF[0] != ZERO)
960                WAXPY(N, rkF[0], Z1, 1, DZ4, 1);
961           if (rkF[1] != ZERO)
962                WAXPY(N, rkF[1], Z2, 1, DZ4, 1);
963           if (rkF[2] != ZERO)
964                WAXPY(N, rkF[2], Z3, 1, DZ4, 1);
965           for (i = 0; i < N; i++)
966                 TMP[i] = Y[i] + Z4[i];
967           FUN_CHEM(*T+H, TMP, DZ1);
968           WAXPY(N, H*rkBgam[4], DZ1, 1, DZ4, 1);
969        }
970        else {
971           /* DZ4(1,N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4; */
972           if (rkD[0] != ZERO)
973                WAXPY(N, rkD[0], Z1, 1, DZ4, 1);
974           if (rkD[1] != ZERO)
975                WAXPY(N, rkD[1], Z2, 1, DZ4, 1);
976           if (rkD[2] != ZERO)
977                WAXPY(N, rkD[2], Z3, 1, DZ4, 1);
978           WAXPY(N, -ONE, Z4, 1, DZ4, 1);
979        } /* end if */
980        Err = RK_ErrorNorm(N,SCAL,DZ4);
981      }
982      else
983      {
984        RK_ErrorEstimate(N,H,*T,Y,FO,E1,IP1,Z1,Z2,Z3,SCAL,&Err,
985                        FirstStep,Reject,ISTATUS);
986      } /* end if */
987
988      /*~~~> Computation of new step size Hnew */
989      Fac = pow(Err, (-ONE/rkELO))
990            *MIN(FacSafe,(ONE+2*NewtonMaxit)/(NewtonIter+2*NewtonMaxit));
991      Fac = MIN(FacMax,MAX(FacMin,Fac));
992      Hnew = Fac*H;
993
994 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
995 /*~~~> Accept/reject step */
996 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
997      // accept:
998      if (Err < ONE) { /*~~~> STEP IS ACCEPTED */
999        FirstStep = 0;
1000        ISTATUS[Nacc]++;
1001        if (Gustafsson == 1) {
1002           /*~~~> Predictive controller of Gustafsson */
1003           if (ISTATUS[Nacc] > 1) {
1004              FacGus = FacSafe*(H/Hacc)*pow(Err*Err/ErrOld,(KPP_REAL)(-0.25));
1005              FacGus = MIN(FacMax,MAX(FacMin,FacGus));
1006              Fac = MIN(Fac,FacGus);   
1007              Hnew = Fac*H;
1008           } /* end if */
1009           Hacc = H;
1010           ErrOld = MAX((KPP_REAL)1.0e-02,Err);
1011        } /* end if */
1012        Hold = H;
1013        *T = *T + H;
1014        /* Update solution: Y <- Y + sum(d_i Z-i) */
1015        if (rkD[0] != ZERO)
1016           WAXPY(N, rkD[0], Z1, 1, Y, 1);
1017        if (rkD[1] != ZERO)
1018           WAXPY(N, rkD[1], Z2, 1, Y, 1);
1019        if (rkD[2] != ZERO)
1020           WAXPY(N, rkD[2], Z3, 1, Y, 1);
1021        /* Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 */
1022        if (StartNewton == 1)
1023           RK_Interpolate("make",N,H,Hold,Z1,Z2,Z3,CONT);
1024        RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL);
1025        RSTATUS[Ntexit] = *T;
1026        RSTATUS[Nhnew] = Hnew;
1027        RSTATUS[Nhacc] = H;
1028        Hnew = Tdirection*MIN( MAX(ABS(Hnew),Hmin), Hmax );
1029        if (Reject == 1)
1030           Hnew = Tdirection*MIN(ABS(Hnew),ABS(H));
1031        Reject = 0;
1032        if ((*T+Hnew/Qmin-Tend)*Tdirection >= ZERO)
1033           H = Tend - *T;
1034        else {
1035           Hratio = Hnew/H;
1036           /* Reuse the LU decomposition */
1037           SkipLU = (Theta <= ThetaMin) && (Hratio >= Qmin) && (Hratio <= Qmax);
1038           if ( SkipLU == 0)
1039                H = Hnew;
1040        } /* end if */
1041        /* If convergence is fast enough, do not update Jacobian */
1042        /* SkipJac = (Theta <= ThetaMin); */
1043        SkipJac = 0;
1044      }
1045
1046      /*~~~> Step is rejected */
1047      else {
1048        if ((FirstStep == 1)  || (Reject == 1))
1049           H = FacRej*H;
1050        else
1051           H = Hnew;
1052        Reject  = 1;
1053        SkipJac = 1; /* Skip if rejected - Jac is independent of H */
1054        SkipLU  = 0;
1055        if (ISTATUS[Nacc] >= 1)
1056           ISTATUS[Nrej]++;
1057      } /* end if accept */
1058   } /* while: time Tloop */
1059 
1060   /*~~~> Successful exit */
1061   *IERR = 1;
1062
1063} /* RK_Integrator */
1064
1065/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1066        Handles all error messages and returns IERR = error Code
1067~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1068void RK_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H, int* IERR)
1069{
1070   Code = *IERR;
1071   printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
1072   printf("\nForced exit from RungeKutta due to the following error:\n");
1073
1074   switch (Code) {
1075   case -1:
1076      printf("--> Improper value for maximal no of steps");
1077      break;
1078   case -2:
1079      printf("--> Improper value for maximal no of Newton iterations");
1080      break;
1081   case -3:
1082      printf("--> Hmin/Hmax/Hstart must be positive");
1083      break;
1084   case -4:
1085      printf("--> Improper values for FacMin/FacMax/FacSafe/FacRej");
1086      break;
1087   case -5:
1088      printf("--> Improper value for ThetaMin");
1089      break;
1090   case -6:
1091      printf("--> Newton stopping tolerance too small");
1092      break;
1093   case -7:
1094      printf("--> Improper values for Qmin, Qmax");
1095      break;
1096   case -8:
1097      printf("--> Tolerances are too small");
1098      break;
1099   case -9:
1100      printf("--> No of steps exceeds maximum bound");
1101      break;
1102   case -10:
1103      printf("--> Step size too small: (T + 10*H = T) or H < Roundoff");
1104      break;
1105   case -11:
1106      printf("--> Matrix is repeatedly singular");
1107      break;
1108   case -12:
1109      printf("--> Non-convergence of Newton iterations");
1110      break;
1111   case -13:
1112      printf("--> Requested RK method not implemented");
1113      break;
1114   default:
1115      printf("Unknown Error code: %d \n", Code);
1116   } /* end switch */ 
1117
1118   printf("\n     T=%e,  H =%e", T, H);
1119   printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
1120
1121} /* RK_ErrorMsg */
1122
1123/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1124/*      Handles all error messages and returns SCAL */
1125/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1126void RK_ErrorScale(int N, 
1127   /*~~~> Input arguments: */
1128        int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[], KPP_REAL Y[],
1129   /*~~~> Output arguments: */
1130        KPP_REAL SCAL[])
1131{
1132   /*printf("Starting RK_ErrorScale\n");*/
1133   int i;
1134   if (ITOL == 0) {
1135        for (i = 0; i < N; i++) {
1136           SCAL[i] = ONE/(AbsTol[0]+RelTol[0]*ABS(Y[i]));
1137        } /* end for loop*/
1138   }
1139   else {
1140        for (i = 0; i< N; i++) {
1141           SCAL[i] = ONE/(AbsTol[i]+RelTol[i]*ABS(Y[i]));
1142        } /* end for loop */
1143   } /* end if */
1144   /*printf("Ending RK_ErrorScale\n");*/
1145} /* RK_ErrorScale */ 
1146
1147/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1148void RK_Transform(int N, KPP_REAL Tr[][3], KPP_REAL Z1[], KPP_REAL Z2[],
1149                KPP_REAL Z3[], KPP_REAL W1[], KPP_REAL W2[], KPP_REAL W3[])
1150-->     W <-- Tr x Z
1151~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1152{
1153   int i;
1154   KPP_REAL x1,x2,x3;
1155   for (i = 0; i < N; i++)
1156   {
1157        x1 = Z1[i];
1158        x2 = Z2[i];
1159        x3 = Z3[i];
1160        W1[i] = Tr[0][0]*x1 + Tr[0][1]*x2 + Tr[0][2]*x3;
1161        W2[i] = Tr[1][0]*x1 + Tr[1][1]*x2 + Tr[1][2]*x3;
1162        W1[i] = Tr[2][0]*x1 + Tr[2][1]*x2 + Tr[2][2]*x3;
1163   }  end for loop
1164}  end RK_Transform */
1165
1166/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1167/*--> Constructs or evaluates a quadratic polynomial that
1168      interpolates the Z solution at current step and
1169      provides starting values for the next step
1170~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1171void RK_Interpolate(char action[], int N, KPP_REAL H, KPP_REAL Hold, KPP_REAL Z1[],
1172        KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL CONT[][RKmax])
1173{
1174   int i;
1175   KPP_REAL r,x1,x2,x3,den;
1176 
1177   /* Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 */
1178   if (action == "make") {
1179        den = (rkC[2]-rkC[1])*(rkC[1]-rkC[0])*(rkC[0]-rkC[2]);
1180        for (i = 0; i < N; i++)
1181        {
1182           CONT[i][0] = (-rkC[2]*rkC[2]*rkC[1]*Z1[i]
1183                        +Z3[i]*rkC[1]*rkC[0]*rkC[0]
1184                        +rkC[1]*rkC[1]*rkC[2]*Z1[i]
1185                        -rkC[1]*rkC[1]*rkC[0]*Z3[i]
1186                        +rkC[2]*rkC[2]*rkC[0]*Z2[i]
1187                        -Z2[i]*rkC[2]*rkC[0]*rkC[0])/den-Z3[i];
1188           CONT[i][1] = -( rkC[0]*rkC[0]*(Z3[i]-Z2[i])
1189                        + rkC[1]*rkC[1]*(Z1[i]-Z3[i])
1190                        + rkC[2]*rkC[2]*(Z2[i]-Z1[i]) )/den;
1191           CONT[i][2] = ( rkC[0]*(Z3[i]-Z2[i])
1192                        + rkC[1]*(Z1[i]-Z3[i])
1193                        + rkC[2]*(Z2[i]-Z1[i]) )/den;
1194        } /* end for loop */
1195   }
1196   /* Evaluate quadratic polynomial */
1197   else if (action == "eval") {
1198        r  = H / Hold;
1199        x1 = ONE + rkC[0]*r;
1200        x2 = ONE + rkC[1]*r;
1201        x3 = ONE + rkC[2]*r;
1202        for (i = 0; i < N; i++)
1203        {
1204           Z1[i] = CONT[i][0]+x1*(CONT[i][1]+x1*CONT[i][2]);   
1205           Z2[i] = CONT[i][0]+x2*(CONT[i][1]+x2*CONT[i][2]);   
1206           Z3[i] = CONT[i][0]+x3*(CONT[i][1]+x3*CONT[i][2]);
1207        } /* end for loop */
1208   } /* end if */
1209} /* RK_Interpolate */
1210
1211/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1212/*~~~> Prepare the right-hand side for Newton iterations
1213       R = Z - hA x F
1214~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1215void RK_PrepareRHS(int N, KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL FO[],
1216        KPP_REAL Z1[], KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL R1[], KPP_REAL R2[],
1217        KPP_REAL R3[])
1218{
1219   KPP_REAL TMP[N], F[N];
1220   WCOPY(N,Z1,1,R1,1); /* R1 <- Z1 */   
1221   WCOPY(N,Z2,1,R2,1); /* R2 <- Z2 */   
1222   WCOPY(N,Z3,1,R3,1); /* R3 <- Z3 */   
1223 
1224   if (rkMethod == L3A)
1225   {
1226        WAXPY(N,-H*rkA[0][0],FO,1,R1,1); /* R1 <- R1 - h*A_10*FO */
1227        WAXPY(N,-H*rkA[1][0],FO,1,R2,1); /* R2 <- R2 - h*A_20*FO */
1228        WAXPY(N,-H*rkA[2][0],FO,1,R3,1); /* R3 <- R3 - h*A_30*FO */
1229   } /* end if */
1230
1231   WADD(N,Y,Z1,TMP);                    /* TMP <- Y + Z1 */
1232   FUN_CHEM(T+rkC[0]*H,TMP,F);          /* F1 <- Fun(Y+Z1) */
1233   WAXPY(N,-H*rkA[0][0],F,1,R1,1);      /* R1 <- R1 - h*A_11*F1 */
1234   WAXPY(N,-H*rkA[1][0],F,1,R2,1);      /* R2 <- R2 - h*A_21*F1 */
1235   WAXPY(N,-H*rkA[2][0],F,1,R3,1);      /* R3 <- R3 - h*A_31*F1 */
1236
1237   WADD(N,Y,Z2,TMP);                    /* TMP <- Y + Z2 */
1238   FUN_CHEM(T+rkC[1]*H,TMP,F);          /* F2 <- Fun(Y+Z2) */
1239   WAXPY(N,-H*rkA[0][1],F,1,R1,1);      /* R1 <- R1 - h*A_12*F2 */
1240   WAXPY(N,-H*rkA[1][1],F,1,R2,1);      /* R2 <- R2 - h*A_22*F2 */
1241   WAXPY(N,-H*rkA[2][1],F,1,R3,1);      /* R3 <- R3 - h*A_32*F2 */
1242
1243   WADD(N,Y,Z3,TMP);                    /* TMP <- Y + Z3 */
1244   FUN_CHEM(T+rkC[2]*H,TMP,F);          /* F3 <- Fun(Y+Z3) */
1245   WAXPY(N,-H*rkA[0][2],F,1,R1,1);      /* R1 <- R1 - h*A_13*F3 */
1246   WAXPY(N,-H*rkA[1][2],F,1,R2,1);      /* R2 <- R2 - h*A_23*F3 */
1247   WAXPY(N,-H*rkA[2][2],F,1,R3,1);      /* R3 <- R3 - h*A_33*F3 */
1248} /* RK_PrepareRHS */
1249
1250/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1251/*~~~> Compute the matrices E1 and E2 and their decompositions
1252~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1253void RK_Decomp(int N, KPP_REAL H, KPP_REAL FJAC[], KPP_REAL E1[],
1254        int IP1[], KPP_REAL E2R[], KPP_REAL E2I[], int IP2[], 
1255        int* ISING, int ISTATUS[]) 
1256{
1257   /*printf("Starting RK_Decomp\n");*/
1258   KPP_REAL Alpha, Beta, Gamma;
1259   int i,j;
1260
1261   Gamma = rkGamma / H;
1262   Alpha = rkAlpha / H;
1263   Beta  = rkBeta / H;
1264
1265   for (i = 0; i < LU_NONZERO; i++) {
1266        E1[i] = -FJAC[i];
1267   }
1268   for (i = 0; i < NVAR; i++) {
1269        j = LU_DIAG[i];
1270        E1[j] = E1[j] + Gamma;
1271   }
1272   *ISING = KppDecomp(E1);
1273   /*~~~> Note: for a full matrix use Lapack:
1274   for (j = 0; j < N; j++) {
1275        for (i = 0; i < N; i++)
1276        { 
1277           E1[i,j] = -FJAC[i][j];
1278        }
1279        E1[i][j] = E1[i][j]+Gamma;
1280   }
1281   DGETRF(N, N, E1, N, IP1, ISING); */
1282
1283   if (*ISING != 0) {
1284        /*printf("Matrix is singular, ISING=%d\n", *ISING);*/
1285        ISTATUS[Ndec]++;
1286        return;
1287   }
1288
1289   for (i = 0; i < LU_NONZERO; i++) {
1290        E2R[i] = (KPP_REAL)(-FJAC[i]);
1291        E2I[i] = ZERO;
1292   }
1293   for (i = 0; i < NVAR; i++) {
1294        j = LU_DIAG[i];
1295        E2R[j] = E2R[j] + Alpha;
1296        E2I[j] = E2I[j] + Beta;
1297   }
1298   *ISING = KppDecompCmplxR(E2R, E2I);
1299   /*printf("Matrix is singular, ISING=%d\n", *ISING);*/
1300   /*~~~> Note: for a full matrix use Lapack:
1301   for (j = 0; j < N; j++) {
1302        for (i = 0; i < N; i++) {
1303           E2[i,j] = DCMPLX( -FJAC[i,j], ZERO );
1304        }
1305        E2[j,j] = E2[j,j] + CMPLX( Alpha, Beta );
1306   }
1307   ZGETRF(N,N,E2,N,IP2,ISING); */
1308   ISTATUS[Ndec]++;
1309} /*RK_Decomp */
1310
1311/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ 
1312void RK_Solve(int N, KPP_REAL H, KPP_REAL E1[], int IP1[], KPP_REAL E2R[],
1313        KPP_REAL E2I[], int IP2[],KPP_REAL R1[], KPP_REAL R2[], KPP_REAL R3[],
1314        int ISTATUS[])
1315{
1316   KPP_REAL x1, x2, x3;
1317   KPP_REAL BCR[N], BCI[N];
1318   int i;
1319
1320   /* Z <- h^{-1) T^{-1) A^{-1) x Z */
1321   for (i = 0; i < N; i++)
1322   {
1323        x1 = R1[i]/H;
1324        x2 = R2[i]/H;
1325        x3 = R3[i]/H;
1326        R1[i] = rkTinvAinv[0][0]*x1 + rkTinvAinv[0][1]*x2 + rkTinvAinv[0][2]*x3;
1327        R2[i] = rkTinvAinv[1][0]*x1 + rkTinvAinv[1][1]*x2 + rkTinvAinv[1][2]*x3;
1328        R3[i] = rkTinvAinv[2][0]*x1 + rkTinvAinv[2][1]*x2 + rkTinvAinv[2][2]*x3;
1329   } 
1330   KppSolve(E1,R1);
1331   /*~~~> Note: for a full matrix use Lapack:
1332   DGETRS('N',5,1,E1,N,IP1,R1,5,0); */
1333   for(i = 0; i < N; i++)
1334   {
1335        BCR[i] = (KPP_REAL)(R2[i]);
1336        BCI[i] = (KPP_REAL)(R3[i]);
1337   }
1338   KppSolveCmplxR(E2R,E2I,BCR,BCI);
1339   /*~~~> Note: for a full matrix use Lapack:
1340   ZGETRS ('N',N,1,E2,N,IP2,BC,N,0); */
1341   for (i = 0; i < N; i++)
1342   {
1343        R2[i] = (KPP_REAL)( BCR[i] );
1344        R3[i] = (KPP_REAL)( BCI[i] );
1345   }
1346
1347   /* Z <- T x Z */
1348   for (i = 0; i < N; i++)
1349   {
1350        x1 = R1[i];
1351        x2 = R2[i];
1352        x3 = R3[i];
1353        R1[i] = rkT[0][0]*x1 + rkT[0][1]*x2 + rkT[0][2]*x3;
1354        R2[i] = rkT[1][0]*x1 + rkT[1][1]*x2 + rkT[1][2]*x3;
1355        R3[i] = rkT[2][0]*x1 + rkT[2][1]*x2 + rkT[2][2]*x3;
1356   } 
1357   ISTATUS[Nsol]++;
1358} /* RK_Solve */
1359
1360/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1361void RK_ErrorEstimate(int N, KPP_REAL H, KPP_REAL T, KPP_REAL Y[], 
1362        KPP_REAL FO[], KPP_REAL E1[], int IP1[], KPP_REAL Z1[], 
1363        KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL SCAL[], KPP_REAL* Err,
1364        int FirstStep, int Reject, int ISTATUS[])
1365/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1366{
1367   KPP_REAL F1[N],F2[N],TMP[N];
1368   KPP_REAL HrkE1,HrkE2,HrkE3;
1369   int i;
1370
1371   HrkE1  = rkE[0]/H;
1372   HrkE2  = rkE[1]/H;
1373   HrkE3  = rkE[2]/H;
1374
1375   for (i = 0; i < N; i++)
1376   {
1377        F2[i]  = HrkE1*Z1[i]+HrkE2*Z2[i]+HrkE3*Z3[i];
1378        TMP[i] = rkE[0]*FO[i] + F2[i];
1379   }
1380   
1381   KppSolve(E1, TMP);
1382   if ((rkMethod == R1A) || (rkMethod == GAU) || (rkMethod == L3A))
1383        KppSolve(E1,TMP);
1384   if (rkMethod == GAU)
1385        KppSolve(E1,TMP);
1386   /*~~~> Note: for a full matrix use Lapack:
1387        DGETRS ('N',N,1,E1,N,IP1,TMP,N,0);
1388        if ((rkMethod == R1A) || (rkMethod == GAU) || (rkMethod == L3A))
1389           DGETRS('N',N,1,E1,N,IP1,TMP,N,0);
1390        if (rkMethod == GAU)
1391           DGETRS('N',N,1,E1,N,IP1,TMP,N,0); */
1392
1393   *Err = RK_ErrorNorm(N,SCAL,TMP);
1394
1395   if (*Err < ONE)
1396        return;
1397   /* firej */
1398   if ((FirstStep == 1) || (Reject == 1)) {
1399        for (i = 0; i < N; i++) {
1400           TMP[i] = Y[i] + TMP[i];
1401        } /* end for */
1402        FUN_CHEM(T,TMP,F1);
1403        ISTATUS[Nfun]++;
1404        for (i = 0; i < N; i++) {
1405           TMP[i] = F1[i] + F2[i];
1406        } /* end for */
1407
1408        KppSolve(E1, TMP);
1409        /*~~~> Note: for a full matrix use Lapack:
1410           DGETRS ('N',N,1,E1,N,IP1,TMP,N,0); */
1411        *Err = RK_ErrorNorm(N,SCAL,TMP);
1412   } /* end if firej */
1413} /* RK_ErrorEstimate */
1414
1415/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1416   KPP_REAL RK_ErrorNorm(int N, KPP_REAL SCAL[], KPP_REAL DY[])
1417/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1418{
1419   int i;
1420   KPP_REAL RK_ErrorNorm = ZERO;
1421   for (i = 0; i < N; i++) {
1422        RK_ErrorNorm = RK_ErrorNorm + (DY[i]*SCAL[i]) * (DY[i]*SCAL[i]);
1423   }   
1424   RK_ErrorNorm = MAX( SQRT(RK_ErrorNorm/N), (KPP_REAL)1.0e-10 );
1425   return RK_ErrorNorm;
1426} /* RK_ErrorNorm */
1427
1428/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1429        The coefficients of the 3-stage Radau-2A method
1430        (given to ~30 accurate digits)
1431~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1432void Radau2A_Coefficients()
1433{
1434   /*printf("Starting Radau2A\n");*/
1435   /* The coefficients of the Radau2A method */
1436   KPP_REAL b0;
1437   /* b0 = (KPP_REAL)1.0; */
1438   if (SdirkError == 1) {
1439        b0 = (KPP_REAL)(0.2e-01);
1440        /*printf("b0 = %f\n", b0);*/
1441   }
1442   else {
1443        b0 = (KPP_REAL)(0.5e-01);
1444        /*printf("b0 = %f\n", b0);*/
1445   } /*end if */
1446   /* The coefficients of the Radau2A method */
1447   rkMethod = R2A;
1448
1449   rkA[0][0] = (KPP_REAL)(1.968154772236604258683861429918299e-01);
1450   rkA[0][1] = (KPP_REAL)(-6.55354258501983881085227825696087e-02);
1451   rkA[0][2] = (KPP_REAL)(2.377097434822015242040823210718965e-02);
1452   rkA[1][0] = (KPP_REAL)(3.944243147390872769974116714584975e-01);
1453   rkA[1][1] = (KPP_REAL)(2.920734116652284630205027458970589e-01);
1454   rkA[1][2] = (KPP_REAL)(-4.154875212599793019818600988496743e-02);
1455   rkA[2][0] = (KPP_REAL)(3.764030627004672750500754423692808e-01);
1456   rkA[2][1] = (KPP_REAL)(5.124858261884216138388134465196080e-01);
1457   rkA[2][2] = (KPP_REAL)(1.111111111111111111111111111111111e-01);
1458
1459   rkB[0] = (KPP_REAL)(3.764030627004672750500754423692808e-01);
1460   rkB[1] = (KPP_REAL)(5.124858261884216138388134465196080e-01);
1461   rkB[2] = (KPP_REAL)(1111111111111111111111111111111111e-01);
1462
1463   rkC[0] = (KPP_REAL)(1.550510257216821901802715925294109e-01);
1464   rkC[1] = (KPP_REAL)(6.449489742783178098197284074705891e-01);
1465   rkC[2] = ONE;
1466
1467   /* New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j */
1468   rkD[0] = ZERO;
1469   rkD[1] = ZERO;
1470   rkD[2] = ONE;
1471   
1472   /* Classical error estimator: */
1473   /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1474   rkE[0] = ONE*b0;
1475   rkE[1] = (KPP_REAL)(-10.04880939982741556246032950764708)*b0;
1476   rkE[2] = (KPP_REAL)1.382142733160748895793662840980412*b0;
1477   rkE[3] = (KPP_REAL)(-.3333333333333333333333333333333333)*b0;
1478
1479   /* Sdirk error estimator */
1480   rkBgam[0] = b0;
1481   rkBgam[1] = (KPP_REAL).3764030627004672750500754423692807
1482                -(KPP_REAL)1.558078204724922382431975370686279*b0;
1483   rkBgam[2] = (KPP_REAL).8914115380582557157653087040196118*b0
1484                +(KPP_REAL).5124858261884216138388134465196077;
1485   rkBgam[3] = (KPP_REAL)(-.1637777184845662566367174924883037)
1486                -(KPP_REAL).3333333333333333333333333333333333*b0;
1487   rkBgam[4] = (KPP_REAL).2748888295956773677478286035994148;
1488
1489   /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1490   rkTheta[0] = (KPP_REAL)(-1.520677486405081647234271944611547)
1491                -(KPP_REAL)10.04880939982741556246032950764708*b0;
1492   rkTheta[1] = (KPP_REAL)2.070455145596436382729929151810376
1493                +(KPP_REAL)1.382142733160748895793662840980413*b0;
1494   rkTheta[2] = (KPP_REAL)(-.3333333333333333333333333333333333)*b0
1495                -(KPP_REAL).3744441479783868387391430179970741;
1496
1497   /* Local order of error estimator */
1498   if ( b0 == ZERO)
1499        rkELO  = (KPP_REAL)6.0;
1500   else
1501        rkELO  = (KPP_REAL)4.0;
1502
1503   /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1504    ~~~> Diagonalize the RK matrix:
1505         rkTinv * inv(rkA) * rkT =
1506                   |  rkGamma      0           0     |
1507                   |      0      rkAlpha   -rkBeta   |
1508                   |      0      rkBeta     rkAlpha  |
1509     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1510
1511   rkGamma = (KPP_REAL)3.637834252744495732208418513577775;
1512   rkAlpha = (KPP_REAL)2.681082873627752133895790743211112;
1513   rkBeta  = (KPP_REAL)3.050430199247410569426377624787569;
1514
1515   rkT[0][0] = (KPP_REAL)(9.443876248897524148749007950641664e-02);
1516   rkT[0][1] = (KPP_REAL)(-1.412552950209542084279903838077973e-01);
1517   rkT[0][2] = (KPP_REAL)(-3.00291941051474244918611170890539e-02);
1518   rkT[1][0] = (KPP_REAL)(2.502131229653333113765090675125018e-01);
1519   rkT[1][1] = (KPP_REAL)(2.041293522937999319959908102983381e-01);
1520   rkT[1][2] = (KPP_REAL)(3.829421127572619377954382335998733e-01);
1521   rkT[2][0] =  ONE;
1522   rkT[2][1] =  ONE;
1523   rkT[2][2] =  ZERO;
1524
1525   rkTinv[0][0] = (KPP_REAL)4.178718591551904727346462658512057;
1526   rkTinv[0][1] = (KPP_REAL)(3.27682820761062387082533272429617e-01);
1527   rkTinv[0][2] = (KPP_REAL)(5.233764454994495480399309159089876e-01);
1528   rkTinv[1][0] = (KPP_REAL)(-4.178718591551904727346462658512057);
1529   rkTinv[1][1] = (KPP_REAL)(-3.27682820761062387082533272429617e-01);
1530   rkTinv[1][2] = (KPP_REAL)(4.766235545005504519600690840910124e-01);
1531   rkTinv[2][0] = (KPP_REAL)(-5.02872634945786875951247343139544e-01);
1532   rkTinv[2][1] = (KPP_REAL)2.571926949855605429186785353601676;
1533   rkTinv[2][2] = (KPP_REAL)(-5.960392048282249249688219110993024e-01);
1534
1535   rkTinvAinv[0][0] = (KPP_REAL)(1.520148562492775501049204957366528e+01);
1536   rkTinvAinv[0][1] = (KPP_REAL)1.192055789400527921212348994770778;
1537   rkTinvAinv[0][2] = (KPP_REAL)1.903956760517560343018332287285119;
1538   rkTinvAinv[1][0] = (KPP_REAL)(-9.669512977505946748632625374449567);
1539   rkTinvAinv[1][1] = (KPP_REAL)(-8.724028436822336183071773193986487);
1540   rkTinvAinv[1][2] = (KPP_REAL)3.096043239482439656981667712714881;
1541   rkTinvAinv[2][0] = (KPP_REAL)(-1.409513259499574544876303981551774e+01);
1542   rkTinvAinv[2][1] = (KPP_REAL)5.895975725255405108079130152868952;
1543   rkTinvAinv[2][2] = (KPP_REAL)(-1.441236197545344702389881889085515e-01);
1544
1545   rkAinvT[0][0] = (KPP_REAL).3435525649691961614912493915818282;
1546   rkAinvT[0][1] = (KPP_REAL)(-.4703191128473198422370558694426832);
1547   rkAinvT[0][2] = (KPP_REAL).3503786597113668965366406634269080;
1548   rkAinvT[1][0] = (KPP_REAL).9102338692094599309122768354288852;
1549   rkAinvT[1][1] = (KPP_REAL)1.715425895757991796035292755937326;
1550   rkAinvT[1][2] = (KPP_REAL).4040171993145015239277111187301784;
1551   rkAinvT[2][0] = (KPP_REAL)3.637834252744495732208418513577775;
1552   rkAinvT[2][1] = (KPP_REAL)2.681082873627752133895790743211112;
1553   rkAinvT[2][2] = (KPP_REAL)(-3.050430199247410569426377624787569);
1554   /*printf("Ending Radau2A");*/
1555} /* Radau2A_Coefficients */
1556
1557/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1558        The coefficients of the 3-stage Lobatto-3C method
1559        (given to ~30 accurate digits)
1560        The parameter b0 can be chosen to tune the error estimator
1561~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1562void Lobatto3C_Coefficients ()
1563{
1564   KPP_REAL b0; 
1565   rkMethod = L3C;
1566   /* b0 = 1.0d0 */
1567   if (SdirkError == 1)
1568        b0 = (KPP_REAL)0.2;
1569   else
1570        b0 = (KPP_REAL)0.5;
1571   
1572   /* The coefficients of the Lobatto3C method */
1573   rkA[0][0] = (KPP_REAL).1666666666666666666666666666666667;
1574   rkA[0][1] = (KPP_REAL)(-.3333333333333333333333333333333333);
1575   rkA[0][2] = (KPP_REAL).1666666666666666666666666666666667;
1576   rkA[1][0] = (KPP_REAL).1666666666666666666666666666666667;
1577   rkA[1][1] = (KPP_REAL).4166666666666666666666666666666667;
1578   rkA[1][2] = (KPP_REAL)(-.8333333333333333333333333333333333e-01);
1579   rkA[2][0] = (KPP_REAL).1666666666666666666666666666666667;
1580   rkA[2][1] = (KPP_REAL).6666666666666666666666666666666667;
1581   rkA[2][2] = (KPP_REAL).1666666666666666666666666666666667;
1582
1583   rkB[0] = (KPP_REAL).1666666666666666666666666666666667;
1584   rkB[1] = (KPP_REAL).6666666666666666666666666666666667;
1585   rkB[2] = (KPP_REAL).1666666666666666666666666666666667;
1586
1587   rkC[0] = ZERO;
1588   rkC[1] = (KPP_REAL)0.5;
1589   rkC[2] = ONE;
1590
1591   /* Classical error estimator, embedded solution: */
1592   rkBhat[0] = b0;
1593   rkBhat[1] = (KPP_REAL).16666666666666666666666666666666667-b0;
1594   rkBhat[2] = (KPP_REAL).66666666666666666666666666666666667;
1595   rkBhat[3] = (KPP_REAL).16666666666666666666666666666666667;
1596
1597   /* New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j */
1598   rkD[0] = ZERO;
1599   rkD[1] = ZERO;
1600   rkD[2] = ONE;
1601
1602   /* Classical error estimator: */
1603   /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1604   rkE[0] = (KPP_REAL).3808338772072650364017425226487022*b0;
1605   rkE[1] = (KPP_REAL)(-1.142501631621795109205227567946107)*b0;
1606   rkE[2] = (KPP_REAL)(-1.523335508829060145606970090594809)*b0;
1607   rkE[3] = (KPP_REAL).3808338772072650364017425226487022*b0;
1608
1609   /* Sdirk error estimator */
1610   rkBgam[0] = b0;
1611   rkBgam[1] = (KPP_REAL).1666666666666666666666666666666667-1.0*b0;
1612   rkBgam[2] = (KPP_REAL).6666666666666666666666666666666667;
1613   rkBgam[3] = (KPP_REAL)(-.2141672105405983697350758559820354);
1614   rkBgam[4] = (KPP_REAL).3808338772072650364017425226487021;
1615
1616   /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1617   rkTheta[0] = -3.0*b0-(KPP_REAL).3808338772072650364017425226487021;
1618   rkTheta[1] = -4.0*b0+(KPP_REAL)1.523335508829060145606970090594808;
1619   rkTheta[2] = (KPP_REAL)(-.142501631621795109205227567946106)+b0;
1620
1621   /* Local order of error estimator */
1622   if (b0 == ZERO)
1623        rkELO  = 5.0;
1624   else
1625        rkELO  = 4.0;
1626
1627   /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1628   ~~~> Diagonalize the RK matrix:
1629    rkTinv * inv(rkA) * rkT =
1630              |  rkGamma      0           0     |
1631              |      0      rkAlpha   -rkBeta   |
1632              |      0      rkBeta     rkAlpha  |
1633   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1634
1635   rkGamma = (KPP_REAL)2.625816818958466716011888933765284;
1636   rkAlpha = (KPP_REAL)1.687091590520766641994055533117359;
1637   rkBeta  = (KPP_REAL)2.508731754924880510838743672432351;
1638
1639   rkT[0][1] = ONE;
1640   rkT[0][1] = ONE;
1641   rkT[0][2] = ZERO;
1642   rkT[1][0] = (KPP_REAL).4554100411010284672111720348287483;
1643   rkT[1][1] = (KPP_REAL)(-.6027050205505142336055860174143743);
1644   rkT[1][2] = (KPP_REAL)(-.4309321229203225731070721341350346);
1645   rkT[2][0] = (KPP_REAL)2.195823345445647152832799205549709;
1646   rkT[2][1] = (KPP_REAL)(-1.097911672722823576416399602774855);
1647   rkT[2][2] = (KPP_REAL).7850032632435902184104551358922130;
1648
1649   rkTinv[0][0] = (KPP_REAL).4205559181381766909344950150991349;
1650   rkTinv[0][1] = (KPP_REAL).3488903392193734304046467270632057;
1651   rkTinv[0][2] = (KPP_REAL).1915253879645878102698098373933487;
1652   rkTinv[1][0] = (KPP_REAL).5794440818618233090655049849008650;
1653   rkTinv[1][1] = (KPP_REAL)(-.3488903392193734304046467270632057);
1654   rkTinv[1][2] = (KPP_REAL)(-.1915253879645878102698098373933487);
1655   rkTinv[2][0] = (KPP_REAL)(-.3659705575742745254721332009249516);
1656   rkTinv[2][1] = (KPP_REAL)(-1.463882230297098101888532803699806);
1657   rkTinv[2][2] = (KPP_REAL).4702733607340189781407813565524989;
1658
1659   rkTinvAinv[0][0] = (KPP_REAL)1.104302803159744452668648155627548;
1660   rkTinvAinv[0][1] = (KPP_REAL).916122120694355522658740710823143;
1661   rkTinvAinv[0][2] = (KPP_REAL).5029105849749601702795812241441172;
1662   rkTinvAinv[1][0] = (KPP_REAL)1.895697196840255547331351844372453;
1663   rkTinvAinv[1][1] = (KPP_REAL)3.083877879305644477341259289176857;
1664   rkTinvAinv[1][2] = (KPP_REAL)(-1.502910584974960170279581224144117);
1665   rkTinvAinv[2][0] = (KPP_REAL).8362439183082935036129145574774502;
1666   rkTinvAinv[2][1] = (KPP_REAL)(-3.344975673233174014451658229909802);
1667   rkTinvAinv[2][2] = (KPP_REAL).312908409479233358005944466882642;
1668
1669   rkAinvT[0][0] = (KPP_REAL)2.625816818958466716011888933765282;
1670   rkAinvT[0][1] = (KPP_REAL)1.687091590520766641994055533117358;
1671   rkAinvT[0][2] = (KPP_REAL)(-2.508731754924880510838743672432351);
1672   rkAinvT[1][0] = (KPP_REAL)1.195823345445647152832799205549710;
1673   rkAinvT[1][1] = (KPP_REAL)(-2.097911672722823576416399602774855);
1674   rkAinvT[1][2] = (KPP_REAL).7850032632435902184104551358922130;
1675   rkAinvT[2][0] = (KPP_REAL)5.765829871932827589653709477334136;
1676   rkAinvT[2][1] = (KPP_REAL).1170850640335862051731452613329320;
1677   rkAinvT[2][2] = (KPP_REAL)4.078738281412060947659653944216779;
1678} /* Lobatto3C_Coefficients */
1679
1680/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1681        The coefficients of the 3-stage Gauss method
1682        (given to ~30 accurate digits)
1683        The parameter b3 can be chosen by the user
1684        to tune the error estimator
1685  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1686void Gauss_Coefficients()
1687{
1688   rkMethod = GAU;
1689   
1690   /* The coefficients of the Gauss method */
1691   KPP_REAL b0;
1692   /*b0 = (KPP_REAL)4.0; */
1693   b0 = (KPP_REAL)0.1;
1694
1695   /* The coefficients of the Gauss method */
1696   rkA[0][0] = (KPP_REAL).1388888888888888888888888888888889;
1697   rkA[0][1] = (KPP_REAL)(-.359766675249389034563954710966045e-01);
1698   rkA[0][2] = (KPP_REAL)(.97894440153083260495800422294756e-02);
1699   rkA[1][0] = (KPP_REAL).3002631949808645924380249472131556;
1700   rkA[1][1] = (KPP_REAL).2222222222222222222222222222222222;
1701   rkA[1][2] = (KPP_REAL)(-.224854172030868146602471694353778e-01);
1702   rkA[2][0] = (KPP_REAL).2679883337624694517281977355483022;
1703   rkA[2][1] = (KPP_REAL).4804211119693833479008399155410489;
1704   rkA[2][2] = (KPP_REAL).1388888888888888888888888888888889;
1705
1706   rkB[0] = (KPP_REAL).2777777777777777777777777777777778;
1707   rkB[1] = (KPP_REAL).4444444444444444444444444444444444;
1708   rkB[2] = (KPP_REAL).2777777777777777777777777777777778;
1709
1710   rkC[0] = (KPP_REAL).1127016653792583114820734600217600;
1711   rkC[1] = (KPP_REAL).5000000000000000000000000000000000;
1712   rkC[2] = (KPP_REAL).8872983346207416885179265399782400;
1713
1714   /* Classical error estimator, embedded solution: */
1715   rkBhat[0] = b0;
1716   rkBhat[1] = (KPP_REAL)(-1.4788305577012361475298775666303999)*b0
1717                 +(KPP_REAL).27777777777777777777777777777777778;
1718   rkBhat[2] = (KPP_REAL).44444444444444444444444444444444444
1719                 +(KPP_REAL).66666666666666666666666666666666667*b0;
1720   rkBhat[3] = (KPP_REAL)(-.18783610896543051913678910003626672)*b0
1721                 +(KPP_REAL).27777777777777777777777777777777778;
1722
1723   /* New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j */
1724   rkD[0] = (KPP_REAL)(.1666666666666666666666666666666667e+01);
1725   rkD[1] = (KPP_REAL)(-.1333333333333333333333333333333333e+01);
1726   rkD[2] = (KPP_REAL)(.1666666666666666666666666666666667e+01);
1727
1728   /* Classical error estimator: */
1729   /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1730   rkE[0] = (KPP_REAL).2153144231161121782447335303806954*b0;
1731   rkE[1] = (KPP_REAL)(-2.825278112319014084275808340593191)*b0;
1732   rkE[2] = (KPP_REAL).2870858974881495709929780405075939*b0;
1733   rkE[3] = (KPP_REAL)(-.4558086256248162565397206448274867e-01)*b0;
1734
1735   /* Sdirk error estimator */
1736   rkBgam[0] = ZERO;
1737   rkBgam[1] = (KPP_REAL).2373339543355109188382583162660537;
1738   rkBgam[2] = (KPP_REAL).5879873931885192299409334646982414;
1739   rkBgam[3] = (KPP_REAL)(-.4063577064014232702392531134499046e-01);
1740   rkBgam[4] = (KPP_REAL).2153144231161121782447335303806955;
1741
1742   /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1743   rkTheta[0] = (KPP_REAL)(-2.594040933093095272574031876464493);
1744   rkTheta[1] = (KPP_REAL)1.824611539036311947589425112250199;
1745   rkTheta[2] = (KPP_REAL).1856563166634371860478043996459493;
1746
1747   /* ELO = local order of classical error estimator */
1748   rkELO = (KPP_REAL)4.0;
1749
1750   /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1751   ~~~> Diagonalize the RK matrix:
1752    rkTinv * inv(rkA) * rkT =
1753              |  rkGamma      0           0     |
1754              |      0      rkAlpha   -rkBeta   |
1755              |      0      rkBeta     rkAlpha  |
1756   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1757
1758   rkGamma = (KPP_REAL)4.644370709252171185822941421408064;
1759   rkAlpha = (KPP_REAL)3.677814645373914407088529289295970;
1760   rkBeta  = (KPP_REAL)3.508761919567443321903661209182446;
1761
1762   rkT[0][0] = (KPP_REAL)(.7215185205520017032081769924397664e-01);
1763   rkT[0][1] = (KPP_REAL)(-.8224123057363067064866206597516454e-01);
1764   rkT[0][2] = (KPP_REAL)(-.6012073861930850173085948921439054e-01);
1765   rkT[1][0] = (KPP_REAL).1188325787412778070708888193730294;
1766   rkT[1][1] = (KPP_REAL)(.5306509074206139504614411373957448e-01);
1767   rkT[1][2] = (KPP_REAL).3162050511322915732224862926182701;
1768   rkT[2][0] = ONE;
1769   rkT[2][1] = ONE;
1770   rkT[2][2] = ZERO;
1771
1772   rkTinv[0][0] = (KPP_REAL)5.991698084937800775649580743981285;
1773   rkTinv[0][1] = (KPP_REAL)1.139214295155735444567002236934009;
1774   rkTinv[0][2] = (KPP_REAL).4323121137838583855696375901180497;
1775   rkTinv[1][0] = (KPP_REAL)(-5.991698084937800775649580743981285);
1776   rkTinv[1][1] = (KPP_REAL)(-1.139214295155735444567002236934009);
1777   rkTinv[1][2] = (KPP_REAL).5676878862161416144303624098819503;
1778   rkTinv[2][0] = (KPP_REAL)(-1.246213273586231410815571640493082);
1779   rkTinv[2][1] = (KPP_REAL)2.925559646192313662599230367054972;
1780   rkTinv[2][2] = (KPP_REAL)(-.2577352012734324923468722836888244);
1781
1782   rkTinvAinv[0][0] = (KPP_REAL)27.82766708436744962047620566703329;
1783   rkTinvAinv[0][1] = (KPP_REAL)5.290933503982655311815946575100597;
1784   rkTinvAinv[0][2] = (KPP_REAL)2.007817718512643701322151051660114;
1785   rkTinvAinv[1][0] = (KPP_REAL)(-17.66368928942422710690385180065675);
1786   rkTinvAinv[1][1] = (KPP_REAL)(-14.45491129892587782538830044147713);
1787   rkTinvAinv[1][2] = (KPP_REAL)2.992182281487356298677848948339886;
1788   rkTinvAinv[2][0] = (KPP_REAL)(-25.60678350282974256072419392007303);
1789   rkTinvAinv[2][1] = (KPP_REAL)6.762434375611708328910623303779923;
1790   rkTinvAinv[2][2] = (KPP_REAL)1.043979339483109825041215970036771;
1791
1792   rkAinvT[0][0] = (KPP_REAL).3350999483034677402618981153470483;
1793   rkAinvT[0][1] = (KPP_REAL)(-.5134173605009692329246186488441294);
1794   rkAinvT[0][2] = (KPP_REAL)(.6745196507033116204327635673208923e-01);
1795   rkAinvT[1][0] = (KPP_REAL).5519025480108928886873752035738885;
1796   rkAinvT[1][1] = (KPP_REAL)1.304651810077110066076640761092008;
1797   rkAinvT[1][2] = (KPP_REAL).9767507983414134987545585703726984;
1798   rkAinvT[2][0] = (KPP_REAL)4.644370709252171185822941421408064;
1799   rkAinvT[2][1] = (KPP_REAL)3.677814645373914407088529289295970;
1800   rkAinvT[2][2] = (KPP_REAL)(-3.508761919567443321903661209182446);
1801} /* Gauss_Coefficients */
1802
1803/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1804        The coefficients of the 3-stage Gauss method
1805        (given to ~30 accurate digits)
1806        The parameter b3 can be chosen by the user
1807        to tune the error estimator
1808~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1809void Radau1A_Coefficients()
1810{
1811   /* The coefficients of the Radau1A method */
1812   KPP_REAL b0;
1813   /*b0 = (KPP_REAL)0.3; */
1814   b0 = (KPP_REAL)0.1;
1815
1816   /* The coefficients of the Radau1A method */
1817   rkMethod = R1A;
1818
1819   rkA[0][0] = (KPP_REAL) .1111111111111111111111111111111111;
1820   rkA[0][1] = (KPP_REAL)(-.1916383190435098943442935597058829);
1821   rkA[0][2] = (KPP_REAL)(.8052720793239878323318244859477174e-01);
1822   rkA[1][0] = (KPP_REAL).1111111111111111111111111111111111;
1823   rkA[1][1] = (KPP_REAL).2920734116652284630205027458970589;
1824   rkA[1][2] = (KPP_REAL)(-.481334970546573839513422644787591e-01);
1825   rkA[2][0] = (KPP_REAL).1111111111111111111111111111111111;
1826   rkA[2][1] = (KPP_REAL).5370223859435462728402311533676479;
1827   rkA[2][2] = (KPP_REAL).1968154772236604258683861429918299;
1828
1829   rkB[0] = (KPP_REAL).1111111111111111111111111111111111;
1830   rkB[1] = (KPP_REAL).5124858261884216138388134465196080;
1831   rkB[2] = (KPP_REAL).3764030627004672750500754423692808;
1832
1833   rkC[0] = ZERO;
1834   rkC[1] = (KPP_REAL).3550510257216821901802715925294109;
1835   rkC[2] = (KPP_REAL).8449489742783178098197284074705891;
1836
1837   /* Classical error estimator, embedded solution: */
1838   rkBhat[0] = b0;
1839   rkBhat[1] = (KPP_REAL).11111111111111111111111111111111111-b0;
1840   rkBhat[2] = (KPP_REAL).51248582618842161383881344651960810;
1841   rkBhat[3] = (KPP_REAL).37640306270046727505007544236928079;
1842
1843   /* New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j */
1844   rkD[0] = (KPP_REAL).3333333333333333333333333333333333;
1845   rkD[1] = (KPP_REAL)(-.8914115380582557157653087040196127);
1846   rkD[2] = (KPP_REAL)(1.558078204724922382431975370686279);
1847
1848   /* Classical error estimator: */
1849   /* H* Sum (b_j-bhat_j) f(Z_j) = H*E(0)*F(0) + Sum E_j Z_j */
1850   rkE[0] = (KPP_REAL).2748888295956773677478286035994148*b0;
1851   rkE[1] = (KPP_REAL)(-1.374444147978386838739143017997074)*b0;
1852   rkE[2] = (KPP_REAL)(-1.335337922441686804550326197041126)*b0;
1853   rkE[3] = (KPP_REAL).235782604058977333559011782643466*b0;
1854
1855   /* Sdirk error estimator */
1856   rkBgam[0] = ZERO;
1857   rkBgam[1] = (KPP_REAL)(.1948150124588532186183490991130616e-01);
1858   rkBgam[2] = (KPP_REAL).7575249005733381398986810981093584;
1859   rkBgam[3] = (KPP_REAL)(-.518952314149008295083446116200793e-01);
1860   rkBgam[4] = (KPP_REAL).2748888295956773677478286035994148;
1861
1862   /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1863   rkTheta[0] = (KPP_REAL)(-1.224370034375505083904362087063351);
1864   rkTheta[1] = (KPP_REAL).9340045331532641409047527962010133;
1865   rkTheta[2] = (KPP_REAL).4656990124352088397561234800640929;
1866
1867   /* ELO = local order of classical error estimator */
1868   rkELO = (KPP_REAL)4.0;
1869
1870   /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1871   ~~~> Diagonalize the RK matrix:
1872    rkTinv * inv(rkA) * rkT =
1873              |  rkGamma      0           0     |
1874              |      0      rkAlpha   -rkBeta   |
1875              |      0      rkBeta     rkAlpha  |
1876   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1877
1878   rkGamma = (KPP_REAL)3.637834252744495732208418513577775;
1879   rkAlpha = (KPP_REAL)2.681082873627752133895790743211112;
1880   rkBeta  = (KPP_REAL)3.050430199247410569426377624787569;
1881
1882   rkT[0][0] = (KPP_REAL).424293819848497965354371036408369;
1883   rkT[0][1] = (KPP_REAL)(-.3235571519651980681202894497035503);
1884   rkT[0][2] = (KPP_REAL)(-.522137786846287839586599927945048);
1885   rkT[1][0] = (KPP_REAL)(.57594609499806128896291585429339e-01);
1886   rkT[1][1] = (KPP_REAL)(.3148663231849760131614374283783e-02);
1887   rkT[1][2] = (KPP_REAL).452429247674359778577728510381731;
1888   rkT[2][0] = ONE;
1889   rkT[2][1] = ONE;
1890   rkT[2][2] = ZERO;
1891
1892   rkTinv[0][0] = (KPP_REAL)1.233523612685027760114769983066164;
1893   rkTinv[0][1] = (KPP_REAL)1.423580134265707095505388133369554;
1894   rkTinv[0][2] = (KPP_REAL).3946330125758354736049045150429624;
1895   rkTinv[1][0] = (KPP_REAL)(-1.233523612685027760114769983066164);
1896   rkTinv[1][1] = (KPP_REAL)(-1.423580134265707095505388133369554);
1897   rkTinv[1][2] = (KPP_REAL).6053669874241645263950954849570376;
1898   rkTinv[2][0] = (KPP_REAL)(-.1484438963257383124456490049673414);
1899   rkTinv[2][1] = (KPP_REAL)2.038974794939896109682070471785315;
1900   rkTinv[2][2] = (KPP_REAL)(-.544501292892686735299355831692542e-01);
1901
1902   rkTinvAinv[0][0] = (KPP_REAL)4.487354449794728738538663081025420;
1903   rkTinvAinv[0][1] = (KPP_REAL)5.178748573958397475446442544234494;
1904   rkTinvAinv[0][2] = (KPP_REAL)1.435609490412123627047824222335563;
1905   rkTinvAinv[1][0] = (KPP_REAL)(-2.854361287939276673073807031221493);
1906   rkTinvAinv[1][1] = (KPP_REAL)(-1.003648660720543859000994063139137e+01);
1907   rkTinvAinv[1][2] = (KPP_REAL)1.789135380979465422050817815017383;
1908   rkTinvAinv[2][0] = (KPP_REAL)(-4.160768067752685525282947313530352);
1909   rkTinvAinv[2][1] = (KPP_REAL)1.124128569859216916690209918405860;
1910   rkTinvAinv[2][2] = (KPP_REAL)1.700644430961823796581896350418417;
1911
1912   rkAinvT[0][0] = (KPP_REAL)1.543510591072668287198054583233180;
1913   rkAinvT[0][1] = (KPP_REAL)(-2.460228411937788329157493833295004);
1914   rkAinvT[0][2] = (KPP_REAL)(-.412906170450356277003910443520499);
1915   rkAinvT[1][0] = (KPP_REAL).209519643211838264029272585946993;
1916   rkAinvT[1][1] = (KPP_REAL)1.388545667194387164417459732995766;
1917   rkAinvT[1][2] = (KPP_REAL)1.20339553005832004974976023130002;
1918   rkAinvT[2][0] = (KPP_REAL)3.637834252744495732208418513577775;
1919   rkAinvT[2][1] = (KPP_REAL)2.681082873627752133895790743211112;
1920   rkAinvT[2][2] = (KPP_REAL)(-3.050430199247410569426377624787569);
1921} /* Radau1A_Coefficients */
1922
1923/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1924        The coefficients of the 4-stage Lobatto-3A method
1925        (given to ~30 accurate digits)
1926  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1927void Lobatto3A_Coefficients()
1928{
1929   /* The coefficients of the Lobatto-3A method */
1930   rkMethod = L3A;
1931
1932   rkA[0][0] = ZERO;
1933   rkA[0][1] = ZERO;
1934   rkA[0][2] = ZERO;
1935   rkA[0][3] = ZERO;
1936   rkA[1][0] = (KPP_REAL).11030056647916491413674311390609397;
1937   rkA[1][1] = (KPP_REAL).1896994335208350858632568860939060;
1938   rkA[1][2] = (KPP_REAL)(-.339073642291438837776604807792215e-01);
1939   rkA[1][3] = (KPP_REAL)(.1030056647916491413674311390609397e-01);
1940   rkA[2][0] = (KPP_REAL)(.73032766854168419196590219427239365e-01);
1941   rkA[2][1] = (KPP_REAL).4505740308958105504443271474458881;
1942   rkA[2][2] = (KPP_REAL).2269672331458315808034097805727606;
1943   rkA[2][3] = (KPP_REAL)(-.2696723314583158080340978057276063e-01);
1944   rkA[3][0] = (KPP_REAL)(.83333333333333333333333333333333333e-01);
1945   rkA[3][1] = (KPP_REAL).4166666666666666666666666666666667;
1946   rkA[3][2] = (KPP_REAL).4166666666666666666666666666666667;
1947   rkA[3][3] = (KPP_REAL)(.8333333333333333333333333333333333e-01);
1948
1949   rkB[0] = (KPP_REAL)(.83333333333333333333333333333333333e-01);
1950   rkB[1] = (KPP_REAL).4166666666666666666666666666666667;
1951   rkB[2] = (KPP_REAL).4166666666666666666666666666666667;
1952   rkB[3] = (KPP_REAL)(.8333333333333333333333333333333333e-01);
1953
1954   rkC[0] = ZERO;
1955   rkC[1] = (KPP_REAL).2763932022500210303590826331268724;
1956   rkC[2] = (KPP_REAL).7236067977499789696409173668731276;
1957   rkC[3] = ONE;
1958
1959   /* New solution: H*Sum B_j*f(Z_j) = Sum D_j*Z_j */
1960   rkD[0] = ZERO;
1961   rkD[1] = ZERO;
1962   rkD[2] = ZERO;
1963   rkD[3] = ONE;
1964
1965   /* Classical error estimator, embedded solution: */
1966   rkBhat[0] = (KPP_REAL)(.90909090909090909090909090909090909e-01);
1967   rkBhat[1] = (KPP_REAL).39972675774621371442114262372173276;
1968   rkBhat[2] = (KPP_REAL).43360657558711961891219070961160058;
1969   rkBhat[3] = (KPP_REAL)(.15151515151515151515151515151515152e-01);
1970
1971   /* Classical error estimator: */
1972   /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1973
1974   rkE[0] = (KPP_REAL)(.1957403846510110711315759367097231e-01);
1975   rkE[1] = (KPP_REAL)(-.1986820345632580910316020806676438);
1976   rkE[2] = (KPP_REAL).1660586371214229125096727578826900;
1977   rkE[3] = (KPP_REAL)(-.9787019232550553556578796835486154e-01);
1978
1979   /* Sdirk error estimator: */
1980   rkF[0] = ZERO;
1981   rkF[1] = (KPP_REAL)(-.66535815876916686607437314126436349);
1982   rkF[2] = (KPP_REAL)1.7419302743497277572980407931678409;
1983   rkF[3] = (KPP_REAL)(-1.2918865386966730694684011822841728);
1984
1985   /* ELO = local order of classical error estimator */
1986   rkELO = (KPP_REAL)4.0;
1987
1988   /* Sdirk error estimator: */
1989   rkBgam[0] = (KPP_REAL)(.2950472755430528877214995073815946e-01);
1990   rkBgam[1] = (KPP_REAL).5370310883226113978352873633882769;
1991   rkBgam[2] = (KPP_REAL).2963022450107219354980459699450564;
1992   rkBgam[3] = (KPP_REAL)(-.7815248400375080035021681445218837e-01);
1993   rkBgam[4] = (KPP_REAL).2153144231161121782447335303806956;
1994
1995   /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1996   rkTheta[0] = (KPP_REAL)0.0;
1997   rkTheta[1] = (KPP_REAL)(-.6653581587691668660743731412643631);
1998   rkTheta[2] = (KPP_REAL)1.741930274349727757298040793167842;
1999   rkTheta[3] = (KPP_REAL)(-.291886538696673069468401182284174);
2000
2001
2002   /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2003   ~~~> Diagonalize the RK matrix:
2004    rkTinv * inv(rkA) * rkT =
2005              |  rkGamma      0           0     |
2006              |      0      rkAlpha   -rkBeta   |
2007              |      0      rkBeta     rkAlpha  |
2008   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2009
2010   rkGamma = (KPP_REAL)4.644370709252171185822941421408063;
2011   rkAlpha = (KPP_REAL)3.677814645373914407088529289295968;
2012   rkBeta  = (KPP_REAL)3.508761919567443321903661209182446;
2013
2014   rkT[0][0] = (KPP_REAL)(.5303036326129938105898786144870856e-01);
2015   rkT[0][1] = (KPP_REAL)(-.7776129960563076320631956091016914e-01);
2016   rkT[0][2] = (KPP_REAL)(.6043307469475508514468017399717112e-02);
2017   rkT[1][0] = (KPP_REAL).2637242522173698467283726114649606;
2018   rkT[1][1] = (KPP_REAL).2193839918662961493126393244533346;
2019   rkT[1][2] = (KPP_REAL).3198765142300936188514264752235344;
2020   rkT[2][0] = ONE;
2021   rkT[2][1] = ZERO;
2022   rkT[2][2] = ZERO;
2023
2024   rkTinv[0][0] = (KPP_REAL)7.695032983257654470769069079238553;
2025   rkTinv[0][1] = (KPP_REAL)(-.1453793830957233720334601186354032);
2026   rkTinv[0][2] = (KPP_REAL).6302696746849084900422461036874826;
2027   rkTinv[1][0] = (KPP_REAL)(-7.695032983257654470769069079238553);
2028   rkTinv[1][1] = (KPP_REAL).1453793830957233720334601186354032;
2029   rkTinv[1][2] = (KPP_REAL).3697303253150915099577538963125174;
2030   rkTinv[2][0] = (KPP_REAL)(-1.066660885401270392058552736086173);
2031   rkTinv[2][1] = (KPP_REAL)3.146358406832537460764521760668932;
2032   rkTinv[2][2] = (KPP_REAL)(-.7732056038202974770406168510664738);
2033
2034   rkTinvAinv[0][0] = (KPP_REAL)35.73858579417120341641749040405149;
2035   rkTinvAinv[0][1] = (KPP_REAL)(-.675195748578927863668368190236025);
2036   rkTinvAinv[0][2] = (KPP_REAL)2.927206016036483646751158874041632;
2037   rkTinvAinv[1][0] = (KPP_REAL)(-24.55824590667225493437162206039511);
2038   rkTinvAinv[1][1] = (KPP_REAL)(-10.50514413892002061837750015342036);
2039   rkTinvAinv[1][2] = (KPP_REAL)4.072793983963516353248841125958369;
2040   rkTinvAinv[2][0] = (KPP_REAL)(-30.92301972744621647251975054630589);
2041   rkTinvAinv[2][1] = (KPP_REAL)12.08182467154052413351908559269928;
2042   rkTinvAinv[2][2] = (KPP_REAL)(-1.546411207640594954081233702132946);
2043
2044   rkAinvT[0][0] = (KPP_REAL).2462926658317812882584158369803835;
2045   rkAinvT[0][1] = (KPP_REAL)(-.2647871194157644619747121197289574);
2046   rkAinvT[0][2] = (KPP_REAL).2950720515900466654896406799284586;
2047   rkAinvT[1][0] = (KPP_REAL)1.224833192317784474576995878738004;
2048   rkAinvT[1][1] = (KPP_REAL)1.929224190340981580557006261869763;
2049   rkAinvT[1][2] = (KPP_REAL).4066803323234419988910915619080306;
2050   rkAinvT[2][0] = (KPP_REAL)4.644370709252171185822941421408064;
2051   rkAinvT[2][1] = (KPP_REAL)3.677814645373914407088529289295968;
2052   rkAinvT[2][2] = (KPP_REAL)(-3.508761919567443321903661209182446);
2053} /* Lobatto3A_Coefficients */
2054
2055/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2056        END SUBROUTINE RungeKutta ! and all its internal procedures
2057  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2058
2059/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2060void FUN_CHEM(KPP_REAL T, KPP_REAL V[], KPP_REAL FCT[])
2061{
2062   KPP_REAL Told;
2063
2064   Told = TIME;
2065   TIME = T;
2066   Update_SUN();
2067   Update_RCONST();
2068   Update_PHOTO();
2069   TIME = Told;
2070
2071   Fun(V, FIX, RCONST, FCT);
2072
2073} /* FUN_CHEM */
2074
2075/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2076void JAC_CHEM (KPP_REAL T, KPP_REAL V[], KPP_REAL JF[])
2077{
2078   KPP_REAL Told;
2079
2080   Told = TIME;
2081   TIME = T;
2082   Update_SUN();
2083   Update_RCONST();
2084   Update_PHOTO();
2085   TIME = Told;
2086   Jac_SP(V, FIX, RCONST, JF);
2087} /* JAC_CHEM */
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |