Home


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

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/rosenbrock.c @ 4143

Last change on this file since 4143 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk
File size: 44.5 KB
Line 
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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
90void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
133int 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
441int 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
623KPP_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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
654void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
677char 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
732int 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
770void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
825void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
889void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
967void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
1039void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
1145void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
1175void 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/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/   
1196void 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
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |