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_adj.c @ 4797

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

Merge of branch palm4u into trunk
File size: 80.2 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
6/* Numerical Constants */
7#define ZERO        (KPP_REAL)0.0
8#define ONE         (KPP_REAL)1.0
9#define HALF        (KPP_REAL)0.5
10#define DeltaMin    (KPP_REAL)1.0e-5
11enum boolean { FALSE=0, TRUE=1 };
12
13/*  Statistics on the work performed by the Rosenbrock method */
14enum statistics { Nfun=1, Njac=2, Nstp=3, Nacc=4, Nrej=5, Ndec=6, Nsol=7, 
15                  Nsng=8, Ntexit=1, Nhexit=2, Nhnew=3 };
16
17/*~~~>  Parameters of the Rosenbrock method, up to 6 stages */
18int ros_S, rosMethod;
19enum ros_Params { RS2=1, RS3=2, RS4=3, RD3=4, RD4=5 };
20KPP_REAL ros_A[15], ros_C[15], ros_M[6], ros_E[6], ros_Alpha[6], ros_Gamma[6], 
21       ros_ELO;
22int ros_NewF[6]; /* Holds Boolean values */
23char ros_Name[12]; /* Length 12 */
24
25/*~~~>  Types of Adjoints Implemented */
26enum adjoint { Adj_none=1, Adj_discrete=2, Adj_continuous=3, 
27               Adj_simple_continuous=4 };
28
29/*~~~>  Checkpoints in memory */
30int bufsize = 200000;
31int stack_ptr; /* last written entry */
32KPP_REAL *chk_H, *chk_T;
33KPP_REAL **chk_Y, **chk_K, **chk_J; /* 2D arrays */
34KPP_REAL **chk_dY, **chk_d2Y; /* 2D arrays */
35
36/* Function Headers */
37void INTEGRATE_ADJ(int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], 
38                   KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], 
39                   KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], 
40                   KPP_REAL RCNTRL_U[], int ISTATUS_U[], KPP_REAL RSTATUS_U[]);
41int RosenbrockADJ( KPP_REAL Y[], int NADJ, KPP_REAL Lambda[][NVAR],
42                   KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL AbsTol[], 
43                   KPP_REAL RelTol[], KPP_REAL AbsTol_adj[][NVAR], 
44                   KPP_REAL RelTol_adj[][NVAR], KPP_REAL RCNTRL[], 
45                   int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[] );
46void ros_AllocateDBuffers( int S, int SaveLU );
47void ros_FreeDBuffers( int SaveLU );
48void ros_AllocateCBuffers();
49void ros_FreeCBuffers();
50void ros_DPush( int S, KPP_REAL T, KPP_REAL H, KPP_REAL Ystage[], 
51                KPP_REAL K[], KPP_REAL E[], int P[], int SaveLU );
52void ros_DPop( int S, KPP_REAL* T, KPP_REAL* H, KPP_REAL* Ystage, 
53               KPP_REAL* K, KPP_REAL* E, int* P, int SaveLU );
54void ros_CPush( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], 
55                KPP_REAL d2Y[] );
56void ros_CPop( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], 
57               KPP_REAL d2Y[] );
58int ros_ErrorMsg( int Code, KPP_REAL T, KPP_REAL H);
59int ros_FwdInt (KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL T, 
60                KPP_REAL AbsTol[], KPP_REAL RelTol[], int AdjointType, 
61                KPP_REAL Hmin, KPP_REAL Hstart, KPP_REAL Hmax, 
62                KPP_REAL Roundoff, int ISTATUS[], int Max_no_steps, 
63                KPP_REAL RSTATUS[], int Autonomous, int VectorTol, 
64                KPP_REAL FacMax, KPP_REAL FacMin, KPP_REAL FacSafe, 
65                KPP_REAL FacRej, int SaveLU);
66int ros_DadjInt ( int NADJ, KPP_REAL Lambda[][NVAR], KPP_REAL Tstart, 
67                  KPP_REAL Tend, KPP_REAL T, int SaveLU, int ISTATUS[], 
68                  KPP_REAL Roundoff, int Autonomous);
69int ros_CadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, KPP_REAL Tend,
70                  KPP_REAL T, KPP_REAL AbsTol_adj[][NVAR], 
71                  KPP_REAL RelTol_adj[][NVAR], KPP_REAL RSTATUS[], 
72                  KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, 
73                  KPP_REAL Roundoff, int Max_no_steps, int Autonomous, 
74                  int VectorTol, KPP_REAL FacMax, KPP_REAL FacMin, 
75                  KPP_REAL FacSafe, KPP_REAL FacRej, int ISTATUS[] );
76int ros_SimpleCadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, 
77                        KPP_REAL Tend, KPP_REAL T, int ISTATUS[], 
78                        int Autonomous, KPP_REAL Roundoff );
79KPP_REAL ros_ErrorNorm ( KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], 
80                       KPP_REAL AbsTol[], KPP_REAL RelTol[], int VectorTol );
81void ros_FunTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[], 
82                             KPP_REAL Fcn0[], KPP_REAL dFdT[], int ISTATUS[] );
83void ros_JacTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[], 
84                             KPP_REAL Jac0[], KPP_REAL dJdT[], int ISTATUS[] );
85int ros_PrepareMatrix ( KPP_REAL H, int Direction, KPP_REAL gam, 
86                         KPP_REAL Jac0[], KPP_REAL Ghimj[], int Pivot[], 
87                         int ISTATUS[] );
88void ros_Decomp( KPP_REAL A[], int Pivot[], int* ising, int ISTATUS[] );
89void ros_Solve( char How, KPP_REAL A[], int Pivot[], KPP_REAL b[], 
90                int ISTATUS[] );
91void ros_cadj_Y( KPP_REAL T, KPP_REAL Y[] );
92void ros_Hermite3( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], 
93                   KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Y[] );
94void ros_Hermite5( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], 
95                   KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Ha[], 
96                   KPP_REAL Hb[], KPP_REAL Y[] );
97void Ros2();
98void Ros3();
99void Ros4();
100void Rodas3();
101void Rodas4();
102void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] );
103void HessTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Hes[] );
104void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Fun [] ); 
105void WSCAL( int N, KPP_REAL Alpha, KPP_REAL X[], int incX );
106void WAXPY( int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], 
107            int incY );
108void WCOPY( int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY );
109KPP_REAL WLAMCH( char C );
110void Update_SUN();
111void Update_RCONST();
112void Fun( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] );
113void Jac_SP( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]);
114void Jac_SP_Vec( KPP_REAL Jac[], KPP_REAL Fcn[], KPP_REAL K[] );
115void JacTR_SP_Vec( KPP_REAL Jac[], KPP_REAL Fcn[], KPP_REAL K[] );
116void HessTR_Vec(KPP_REAL Hess[], KPP_REAL U1[], KPP_REAL U2[], KPP_REAL HTU[]);
117void KppSolve( KPP_REAL A[], KPP_REAL b[] );
118int KppDecomp( KPP_REAL A[] );
119void KppSolveTR( KPP_REAL JVS[], KPP_REAL X[], KPP_REAL XX[] );
120void Hessian( KPP_REAL V[], KPP_REAL F[], KPP_REAL RCT[], KPP_REAL Hess[] );
121
122/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
123void INTEGRATE_ADJ( int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], 
124                    KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], 
125                    KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], 
126                    KPP_REAL RCNTRL_U[], int ISTATUS_U[], 
127                    KPP_REAL RSTATUS_U[] ) {
128/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
129
130/*  Local Variables */
131  KPP_REAL RCNTRL[20], RSTATUS[20];
132  int ICNTRL[20], ISTATUS[20], IERR, i;
133
134  for( i = 0; i < 20; i++ ) {
135    ICNTRL[i]  = 0;
136    RCNTRL[i]  = ZERO;
137    ISTATUS[i] = 0;
138    RSTATUS[i] = ZERO;
139  }
140
141/*~~~> fine-tune the integrator:
142  ICNTRL(1) = 0       ! 0 = non-autonomous, 1 = autonomous
143  ICNTRL(2) = 1       ! 0 = scalar, 1 = vector tolerances
144  RCNTRL(3) = STEPMIN ! starting step
145  ICNTRL(3) = 5       ! choice of the method for forward integration
146  ICNTRL(6) = 1       ! choice of the method for continuous adjoint
147  ICNTRL(7) = 2       ! 1=none, 2=discrete, 3=full continuous,
148                        4=simplified continuous adjoint
149  ICNTRL(8) = 1       ! Save fwd LU factorization: 0=*don't* save, 1=save */
150
151/* if optional parameters are given, and if they are >=0, then they overwrite
152   default settings */
153//  if(ICNTRL_U != NULL) {
154//    for(i=0; i<20; i++)
155//      if(ICNTRL_U[i] > 0)
156//        ICNTRL[i] = ICNTRL_U[i];
157//    } /* end for */
158//  } /* end if */
159
160//  if(RCNTRL_U != NULL) {
161//    for(i=0; i<20; i++)
162//      if(RCNTRL_U[i] > 0)
163//        RCNTRL[i] = RCNTRL_U[i];
164//    } /* end for */
165//  } /* end if */
166
167  IERR = RosenbrockADJ( Y, NADJ, Lambda, TIN, TOUT, ATOL, RTOL, ATOL_adj, 
168                        RTOL_adj, RCNTRL, ICNTRL, RSTATUS, ISTATUS );
169
170  if (IERR < 0)
171    printf( "RosenbrockADJ: Unsucessful step at T=%f (IERR=%d)", TIN, IERR );
172
173  STEPMIN = RSTATUS[Nhexit];
174
175/* if optional parameters are given for output
176         copy to them to return information */
177//  if(ISTATUS_U != NULL)
178//    for(i=0; i<20; i++)
179//      ISTATUS_U[i] = ISTATUS[i];
180//  }
181
182//  if(RSTATUS_U != NULL)
183//    for(i=0; i<20; i++)
184//      RSTATUS_U[i] = RSTATUS[i];
185//  }
186
187} /* End of INTEGRATE_ADJ */
188
189/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
190int RosenbrockADJ( KPP_REAL Y[], int NADJ, KPP_REAL Lambda[][NVAR],
191                   KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL AbsTol[], 
192                   KPP_REAL RelTol[], KPP_REAL AbsTol_adj[][NVAR],
193                   KPP_REAL RelTol_adj[][NVAR], KPP_REAL RCNTRL[],
194                   int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[] ) {
195/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
196
197    ADJ = Adjoint of the Tangent Linear Model of a Rosenbrock Method
198
199    Solves the system y'=F(t,y) using a RosenbrockADJ method defined by:
200
201     G = 1/(H*gamma(1)) - Jac(t0,Y0)
202     T_i = t0 + Alpha(i)*H
203     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
204     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
205         gamma(i)*dF/dT(t0, Y0)
206     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
207
208    For details on RosenbrockADJ methods and their implementation consult:
209      E. Hairer and G. Wanner
210      "Solving ODEs II. Stiff and differential-algebraic problems".
211      Springer series in computational mathematics, Springer-Verlag, 1996. 
212    The codes contained in the book inspired this implementation.
213
214    (C)  Adrian Sandu, August 2004
215    Virginia Polytechnic Institute and State University
216    Contact: sandu@cs.vt.edu
217    Revised by Philipp Miehe and Adrian Sandu, May 2006
218    Translation F90 to C by Paul Eller, April 2007
219    This implementation is part of KPP - the Kinetic PreProcessor
220~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
221
222~~~>   INPUT ARGUMENTS:
223
224-     Y[NVAR]    = vector of initial conditions (at T=Tstart)
225      NADJ       -> dimension of linearized system,
226                   i.e. the number of sensitivity coefficients
227-     Lambda[NVAR][NADJ] -> vector of initial sensitivity conditions
228                            (at T=Tstart)
229-    [Tstart,Tend]  = time range of integration
230     (if Tstart>Tend the integration is performed backwards in time)
231-    RelTol, AbsTol = user precribed accuracy
232- void Fun( T, Y, Ydot ) = ODE function,
233                       returns Ydot = Y' = F(T,Y)
234- void Jac( T, Y, Jcb ) = Jacobian of the ODE function,
235                       returns Jcb = dF/dY
236-    ICNTRL[0:9]    = integer inputs parameters
237-    RCNTRL[0:9]    = real inputs parameters
238~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
239
240~~~>     OUTPUT ARGUMENTS:
241
242-    Y[NVAR]    -> vector of final states (at T->Tend)
243-    Lambda[NVAR][NADJ] -> vector of final sensitivities (at T=Tend)
244-    ICNTRL[9:18]   -> integer output parameters
245-    RCNTRL[9:18]   -> real output parameters
246-    IERR       -> job status upon return
247       - succes (positive value) or failure (negative value) -
248           =  1 : Success
249           = -1 : Improper value for maximal no of steps
250           = -2 : Selected RosenbrockADJ method not implemented
251           = -3 : Hmin/Hmax/Hstart must be positive
252           = -4 : FacMin/FacMax/FacRej must be positive
253           = -5 : Improper tolerance values
254           = -6 : No of steps exceeds maximum bound
255           = -7 : Step size too small
256           = -8 : Matrix is repeatedly singular
257~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
258
259~~~>     INPUT PARAMETERS:
260
261    Note: For input parameters equal to zero the default values of the
262       corresponding variables are used.
263
264    ICNTRL[0]   = 1: F = F(y)   Independent of T (AUTONOMOUS)
265              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
266
267    ICNTRL[1]   = 0: AbsTol, RelTol are NVAR-dimensional vectors
268              = 1:  AbsTol, RelTol are scalars
269
270    ICNTRL[2]  -> selection of a particular Rosenbrock method
271        = 0 :  default method is Rodas3
272        = 1 :  method is  Ros2
273        = 2 :  method is  Ros3
274        = 3 :  method is  Ros4
275        = 4 :  method is  Rodas3
276        = 5:   method is  Rodas4
277
278    ICNTRL[3]  -> maximum number of integration steps
279        For ICNTRL[4]=0) the default value of BUFSIZE is used
280
281    ICNTRL[5]  -> selection of a particular Rosenbrock method for the
282                continuous adjoint integration - for cts adjoint it
283                can be different than the forward method ICNTRL(3)
284         Note 1: to avoid interpolation errors (which can be huge!)
285                   it is recommended to use only ICNTRL[6] = 2 or 4
286         Note 2: the performance of the full continuous adjoint
287                   strongly depends on the forward solution accuracy Abs/RelTol
288
289    ICNTRL[6] -> Type of adjoint algorithm
290         = 0 : default is discrete adjoint ( of method ICNTRL[3] )
291         = 1 : no adjoint
292         = 2 : discrete adjoint ( of method ICNTRL[3] )
293         = 3 : fully adaptive continuous adjoint ( with method ICNTRL[6] )
294         = 4 : simplified continuous adjoint ( with method ICNTRL[6] )
295
296    ICNTRL[7]  -> checkpointing the LU factorization at each step:
297        ICNTRL[7]=0 : do *not* save LU factorization (the default)
298        ICNTRL[7]=1 : save LU factorization
299        Note: if ICNTRL[7]=1 the LU factorization is *not* saved
300
301~~~>  Real input parameters:
302
303    RCNTRL[0]  -> Hmin, lower bound for the integration step size
304          It is strongly recommended to keep Hmin = ZERO
305
306    RCNTRL[1]  -> Hmax, upper bound for the integration step size
307
308    RCNTRL[2]  -> Hstart, starting value for the integration step size
309
310    RCNTRL[3]  -> FacMin, lower bound on step decrease factor (default=0.2)
311
312    RCNTRL[4]  -> FacMax, upper bound on step increase factor (default=6)
313
314    RCNTRL[5]  -> FacRej, step decrease factor after multiple rejections
315            (default=0.1)
316
317    RCNTRL[6]  -> FacSafe, by which the new step is slightly smaller
318         than the predicted value  (default=0.9)
319
320    RCNTRL[7]  -> ThetaMin. If Newton convergence rate smaller
321                  than ThetaMin the Jacobian is not recomputed;
322                  (default=0.001)
323
324~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
325
326~~~>     OUTPUT PARAMETERS:
327
328    Note: each call to RosenbrockADJ adds the corrent no. of fcn calls
329      to previous value of ISTATUS(1), and similar for the other params.
330      Set ISTATUS[0:9] = 0 before call to avoid this accumulation.
331
332    ISTATUS[0] = No. of function calls
333    ISTATUS[1] = No. of jacobian calls
334    ISTATUS[2] = No. of steps
335    ISTATUS[3] = No. of accepted steps
336    ISTATUS[4] = No. of rejected steps (except at the beginning)
337    ISTATUS[5] = No. of LU decompositions
338    ISTATUS[6] = No. of forward/backward substitutions
339    ISTATUS[7] = No. of singular matrix decompositions
340
341    RSTATUS[0]  -> Texit, the time corresponding to the
342                   computed Y upon return
343    RSTATUS[1]  -> Hexit, last accepted step before exit
344    For multiple restarts, use Hexit as Hstart in the following run
345~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
346
347/*~~~>  Local variables */
348  KPP_REAL Roundoff, FacMin, FacMax, FacRej, FacSafe;
349  KPP_REAL Hmin, Hmax, Hstart;
350  KPP_REAL Texit=0.0;
351  int i, UplimTol, Max_no_steps=0, IERR;
352  int AdjointType=0, CadjMethod=0;
353  int Autonomous, VectorTol, SaveLU; /* Holds boolean values */
354
355  stack_ptr = -1;
356
357/*~~~>  Initialize statistics */
358  for(i=0; i<20; i++) {
359    ISTATUS[i] = 0;
360    RSTATUS[i] = ZERO;
361  }
362
363/*~~~>  Autonomous or time dependent ODE. Default is time dependent. */
364  Autonomous = !(ICNTRL[0] == 0);
365
366/*~~~> For Scalar tolerances(ICNTRL[1] != 0) the code uses AbsTol[0] and
367         RelTol[0]
368   For Vector tolerances(ICNTRL[1] == 0) the code uses AbsTol[1:NVAR] and
369     RelTol[1:NVAR] */
370
371  if (ICNTRL[1] == 0) {
372    VectorTol = TRUE;
373    UplimTol  = NVAR;
374  }
375  else { 
376    VectorTol = FALSE;
377    UplimTol  = 1;
378  }
379
380/*~~~>   Initialize the particular Rosenbrock method selected */
381  switch( ICNTRL[2] ) {
382    case 0:
383    case 4:
384      Rodas3();
385      break;
386    case 1:
387      Ros2();
388      break;
389    case 2:
390      Ros3();
391      break;
392    case 3:
393      Ros4();
394      break;
395    case 5:
396      Rodas4();
397      break;
398    default:
399      printf( "Unknown Rosenbrock method: ICNTRL[2]=%d", ICNTRL[2] );
400      return ros_ErrorMsg(-2, Tstart, ZERO);
401  } /* End switch */
402
403/*~~~>   The maximum number of steps admitted */
404  if (ICNTRL[3] == 0)
405    Max_no_steps = bufsize - 1;
406  else if (Max_no_steps > 0)
407    Max_no_steps = ICNTRL[3];
408  else {
409    printf("User-selected max no. of steps: ICNTRL[3]=%d",ICNTRL[3] );
410    return ros_ErrorMsg(-1,Tstart,ZERO);
411  }
412
413/*~~~>The particular Rosenbrock method chosen for integrating the cts adjoint*/
414  if (ICNTRL[5] == 0)
415    CadjMethod = 4;
416  else if ( (ICNTRL[5] >= 1) && (ICNTRL[5] <= 5) )
417    CadjMethod = ICNTRL[5];
418  else {
419    printf( "Unknown CADJ Rosenbrock method: ICNTRL[5]=%d", CadjMethod );
420    return ros_ErrorMsg(-2,Tstart,ZERO);
421  }
422
423/*~~~>  Discrete or continuous adjoint formulation */
424  if ( ICNTRL[6] == 0 )
425    AdjointType = Adj_discrete;
426  else if ( (ICNTRL[6] >= 1) && (ICNTRL[6] <= 4) )
427    AdjointType = ICNTRL[6];
428  else {
429    printf( "User-selected adjoint type: ICNTRL[6]=%d", AdjointType );
430    return ros_ErrorMsg(-9,Tstart,ZERO);
431  }
432
433/*~~~> Save or not the forward LU factorization */
434  SaveLU = (ICNTRL[7] != 0);
435
436/*~~~>  Unit roundoff (1+Roundoff>1)  */
437  Roundoff = WLAMCH('E');
438
439/*~~~>  Lower bound on the step size: (positive value) */
440  if (RCNTRL[0] == ZERO)
441    Hmin = ZERO;
442  else if (RCNTRL[0] > ZERO)
443    Hmin = RCNTRL[0];
444  else {
445    printf( "User-selected Hmin: RCNTRL[0]=%f", RCNTRL[0] );
446    return ros_ErrorMsg(-3,Tstart,ZERO);
447  }
448
449/*~~~>  Upper bound on the step size: (positive value) */
450  if (RCNTRL[1] == ZERO)
451    Hmax = ABS(Tend-Tstart);
452  else if (RCNTRL[1] > ZERO)
453    Hmax = MIN(ABS(RCNTRL[1]),ABS(Tend-Tstart));
454  else {
455    printf( "User-selected Hmax: RCNTRL[1]=%f", RCNTRL[1] );
456    return ros_ErrorMsg(-3,Tstart,ZERO);
457  }
458
459/*~~~>  Starting step size: (positive value) */
460  if (RCNTRL[2] == ZERO) {
461    Hstart = MAX(Hmin,DeltaMin);
462  }
463  else if (RCNTRL[2] > ZERO)
464    Hstart = MIN(ABS(RCNTRL[2]),ABS(Tend-Tstart));
465  else {
466    printf( "User-selected Hstart: RCNTRL[2]=%f", RCNTRL[2] );
467    return ros_ErrorMsg(-3,Tstart,ZERO);
468  }
469
470/*~~~>  Step size can be changed s.t.  FacMin < Hnew/Hold < FacMax */
471  if (RCNTRL[3] == ZERO)
472    FacMin = (KPP_REAL)0.2;
473  else if (RCNTRL[3] > ZERO)
474    FacMin = RCNTRL[3];
475  else {
476    printf( "User-selected FacMin: RCNTRL[3]=%f", RCNTRL[3] );
477    return ros_ErrorMsg(-4,Tstart,ZERO);
478  }
479  if (RCNTRL[4] == ZERO)
480    FacMax = (KPP_REAL)6.0;
481  else if (RCNTRL[4] > ZERO)
482    FacMax = RCNTRL[4];
483  else {
484    printf( "User-selected FacMax: RCNTRL[4]=%f", RCNTRL[4] );
485    return ros_ErrorMsg(-4,Tstart,ZERO);
486  }
487
488/*~~~>   FacRej: Factor to decrease step after 2 succesive rejections */
489  if (RCNTRL[5] == ZERO)
490    FacRej = (KPP_REAL)0.1;
491  else if (RCNTRL[5] > ZERO)
492    FacRej = RCNTRL[5];
493  else {
494    printf( "User-selected FacRej: RCNTRL[5]=%f", RCNTRL[5] );
495    return ros_ErrorMsg(-4,Tstart,ZERO);
496  }
497
498/*~~~>   FacSafe: Safety Factor in the computation of new step size */
499  if (RCNTRL[6] == ZERO)
500    FacSafe = (KPP_REAL)0.9;
501  else if (RCNTRL[6] > ZERO)
502    FacSafe = RCNTRL[6];
503  else {
504    printf( "User-selected FacSafe: RCNTRL[6]=%f", RCNTRL[6] );
505    return ros_ErrorMsg(-4,Tstart,ZERO);
506  }
507
508/*~~~>  Check if tolerances are reasonable */
509  for(i=0; i < UplimTol; i++) {
510    if ( (AbsTol[i] <= ZERO) || (RelTol[i] <= (KPP_REAL)10.0*Roundoff)
511         || (RelTol[i] >= (KPP_REAL)1.0) ) {
512      printf( " AbsTol[%d] = %f", i, AbsTol[i] );
513      printf( " RelTol[%d] = %f", i, RelTol[i] );
514      return ros_ErrorMsg(-5,Tstart,ZERO);
515    }
516  }
517
518/*~~~>  Allocate checkpoint space or open checkpoint files */
519  if (AdjointType == Adj_discrete) {
520    ros_AllocateDBuffers( ros_S, SaveLU );
521  }
522  else if ( (AdjointType == Adj_continuous) || 
523            (AdjointType == Adj_simple_continuous) ) {
524    ros_AllocateCBuffers();
525  }
526
527/*~~~>  CALL Forward Rosenbrock method */
528  IERR = ros_FwdInt(Y, Tstart, Tend, Texit, AbsTol, RelTol, AdjointType, Hmin,
529                    Hstart, Hmax, Roundoff, ISTATUS, Max_no_steps, 
530                    RSTATUS, Autonomous, VectorTol, FacMax, FacMin,
531                    FacSafe, FacRej, SaveLU);
532
533  printf( "\n\nFORWARD STATISTICS\n" );
534  printf( "Step=%d Acc=%d Rej=%d Singular=%d\n\n", Nstp, Nacc, Nrej, Nsng );
535
536/*~~~>  If Forward integration failed return */
537  if (IERR<0)
538    return IERR;
539
540/*~~~>   Initialize the particular Rosenbrock method for continuous adjoint */
541  if ( (AdjointType == Adj_continuous) || 
542       (AdjointType == Adj_simple_continuous) ) {
543    switch (CadjMethod) {
544      case 1:
545        Ros2();
546        break;
547      case 2:
548        Ros3();
549        break;
550      case 3:
551        Ros4();
552        break;
553      case 4:
554        Rodas3();
555        break;
556      case 5:
557        Rodas4();
558        break;
559      default:
560        printf( "Unknown Rosenbrock method: ICNTRL[2]=%d", ICNTRL[2] );
561        return ros_ErrorMsg(-2,Tstart,ZERO);
562    }
563  } /* End switch */
564
565  switch( AdjointType ) {
566    case Adj_discrete:
567      IERR = ros_DadjInt (NADJ, Lambda, Tstart, Tend, Texit, SaveLU, ISTATUS, 
568                          Roundoff, Autonomous );
569      break;
570    case Adj_continuous:
571      IERR = ros_CadjInt (NADJ, Lambda, Tend, Tstart, Texit, AbsTol_adj, 
572                          RelTol_adj, RSTATUS, Hmin, Hmax, Hstart, Roundoff,
573                          Max_no_steps, Autonomous, VectorTol, FacMax, FacMin,
574                          FacSafe, FacRej, ISTATUS);
575      break;
576    case Adj_simple_continuous:
577      IERR = ros_SimpleCadjInt (NADJ, Lambda, Tstart, Tend, Texit, ISTATUS, 
578                                Autonomous, Roundoff);
579  } /* End switch for AdjointType */
580
581  printf( "ADJOINT STATISTICS\n" );
582  printf( "Step=%d Acc=%d Rej=%d Singular=%d\n",Nstp,Nacc,Nrej,Nsng );
583
584/*~~~>  Free checkpoint space or close checkpoint files */
585  if (AdjointType == Adj_discrete)
586    ros_FreeDBuffers( SaveLU );
587  else if ( (AdjointType == Adj_continuous) || 
588            (AdjointType == Adj_simple_continuous) )
589    ros_FreeCBuffers();
590
591  return IERR;
592}
593
594/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
595void ros_AllocateDBuffers( int S, int SaveLU ) {
596/*~~~>  Allocate buffer space for discrete adjoint
597~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
598   
599  int i;
600   
601  chk_H = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL));
602  if (chk_H == NULL) {
603    printf("Failed allocation of buffer H");
604    exit(0);
605  }
606   
607  chk_T = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL));
608  if (chk_T == NULL) {
609    printf("Failed allocation of buffer T");
610    exit(0);
611  }
612 
613  chk_Y = (KPP_REAL**) malloc(bufsize * sizeof(KPP_REAL*));
614  if (chk_Y == NULL) {
615    printf("Failed allocation of buffer Y");
616    exit(0);
617  }
618  for(i=0; i<bufsize; i++) {
619    chk_Y[i] = (KPP_REAL*) malloc(NVAR * S * sizeof(KPP_REAL));
620    if (chk_Y[i] == NULL) {
621      printf("Failed allocation of buffer Y");
622      exit(0);
623    }
624  }
625   
626  chk_K = (KPP_REAL**) malloc(bufsize * sizeof(KPP_REAL*));
627  if (chk_K == NULL) {
628    printf("Failed allocation of buffer K");
629    exit(0);
630  }
631  for(i=0; i<bufsize; i++) {
632    chk_K[i] = (KPP_REAL*) malloc(NVAR * S * sizeof(KPP_REAL));
633    if (chk_K == NULL) {
634      printf("Failed allocation of buffer K");
635      exit(0);
636    }
637  }
638 
639  if (SaveLU) {
640    chk_J = (KPP_REAL**) malloc(bufsize * sizeof(KPP_REAL*));
641    if (chk_J == NULL) {
642      printf( "Failed allocation of buffer J");
643      exit(0);
644    }
645#ifdef FULL_ALGEBRA
646    for(i=0; i<bufsize; i++) {
647      chk_J[i] = (KPP_REAL*) malloc(NVAR * NVAR * sizeof(KPP_REAL));
648      if (chk_J == NULL) {
649        printf( "Failed allocation of buffer J");
650        exit(0);
651      }
652    }
653#else
654    for(i=0; i<bufsize; i++) {
655      chk_J[i] = (KPP_REAL*) malloc(LU_NONZERO * sizeof(KPP_REAL));
656      if (chk_J == NULL) {
657        printf( "Failed allocation of buffer J");
658        exit(0);
659      }
660    }
661#endif
662  } 
663
664} /* End of ros_AllocateDBuffers */
665
666/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
667void ros_FreeDBuffers( int SaveLU ) {
668/*~~~>  Deallocate buffer space for discrete adjoint
669~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
670
671  int i;
672
673  free(chk_H);
674  free(chk_T);
675
676  for(i=0; i<bufsize; i++) {
677    free(chk_Y[i]);
678    free(chk_K[i]);
679  }
680  free(chk_Y);
681  free(chk_K);
682
683  if (SaveLU) {
684    for(i=0; i<bufsize; i++)
685      free(chk_J[i]);
686    free(chk_J);
687  }
688
689} /* End of ros_FreeDBuffers */
690
691/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
692void ros_AllocateCBuffers() {
693/*~~~>  Allocate buffer space for continuous adjoint
694~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
695
696  int i;
697
698  chk_H = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL));
699  if (chk_H == NULL) {
700    printf( "Failed allocation of buffer H");
701    exit(0);
702  }
703
704  chk_T = (KPP_REAL*) malloc(bufsize * sizeof(KPP_REAL));
705  if (chk_T == NULL) {
706    printf( "Failed allocation of buffer T");
707    exit(0);
708  }
709
710  chk_Y = (KPP_REAL**) malloc(sizeof(KPP_REAL*) * bufsize);
711  if (chk_Y == NULL) {
712    printf( "Failed allocation of buffer Y");
713    exit(0);
714  }
715  for(i=0; i<bufsize; i++) {
716    chk_Y[i] = (KPP_REAL*) malloc( sizeof(KPP_REAL)* NVAR );
717    if (chk_Y == NULL) {
718      printf( "Failed allocation of buffer Y");
719      exit(0);
720    }
721  }
722
723  chk_dY = (KPP_REAL**) malloc(sizeof(KPP_REAL*) * bufsize);
724  if (chk_dY == NULL) {
725    printf( "Failed allocation of buffer dY");
726    exit(0);
727  }
728  for(i=0; i<bufsize; i++) {
729    chk_dY[i] = (KPP_REAL*) malloc( sizeof(KPP_REAL) * NVAR);
730    if (chk_dY == NULL) {
731      printf( "Failed allocation of buffer dY");
732      exit(0);
733    }
734  }
735
736  chk_d2Y = (KPP_REAL**) malloc(sizeof(KPP_REAL*) * bufsize);
737  if (chk_d2Y == NULL) {
738    printf( "Failed allocation of buffer d2Y");
739    exit(0);
740  }
741  for(i=0; i<bufsize; i++) {
742    chk_d2Y[i] = (KPP_REAL*) malloc( sizeof(KPP_REAL) * NVAR);
743    if (chk_d2Y == NULL) {
744      printf( "Failed allocation of buffer d2Y");
745      exit(0);
746    }
747  }
748} /* End of ros_AllocateCBuffers */
749
750/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
751void ros_FreeCBuffers() {
752/*~~~>  Dallocate buffer space for continuous adjoint
753~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
754 
755  int i;
756
757  free(chk_H);
758  free(chk_T);
759
760  for(i=0; i<bufsize; i++) {
761    free(chk_Y[i]);
762    free(chk_dY[i]);
763    free(chk_d2Y[i]);
764  }
765
766  free(chk_Y);
767  free(chk_dY);
768  free(chk_d2Y);
769
770} /* End of ros_FreeCBuffers */
771
772/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
773void ros_DPush( int S, KPP_REAL T, KPP_REAL H, KPP_REAL Ystage[], 
774                KPP_REAL K[], KPP_REAL E[], int P[], int SaveLU ) {
775/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
776~~~> Saves the next trajectory snapshot for discrete adjoints
777~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
778
779  int i;
780
781  stack_ptr = stack_ptr + 1;
782  if ( stack_ptr >= bufsize ) {
783    printf( "Push failed: buffer overflow" );
784    exit(0);
785  }
786
787  chk_H[ stack_ptr ] = H;
788  chk_T[ stack_ptr ] = T;
789  for(i=0; i<NVAR*S; i++) {
790    chk_Y[stack_ptr][i] = Ystage[i];
791    chk_K[stack_ptr][i] = K[i];
792  }
793
794  if (SaveLU) {
795#ifdef FULL_ALGEBRA
796    int j;
797    for(j=0; j<NVAR; j++) {
798      for(i=0; i<NVAR; i++)
799        chk_J[stack_ptr][i][j] = E[i][j];
800      chk_P[stack_ptr][j] = P[j];
801    }
802#else
803    for(i=0; i<LU_NONZERO; i++)
804      chk_J[stack_ptr][i] = E[i];
805#endif
806  }
807
808} /* End of ros_DPush */
809
810/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
811void ros_DPop( int S, KPP_REAL* T, KPP_REAL* H, KPP_REAL* Ystage,
812               KPP_REAL* K, KPP_REAL* E, int* P, int SaveLU ) {
813/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
814~~~> Retrieves the next trajectory snapshot for discrete adjoints
815~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
816   
817  int i;
818 
819  if ( stack_ptr < 0 ) {
820    printf( "Pop failed: empty buffer" );
821    exit(0);
822  }
823
824  *H = chk_H[ stack_ptr ];
825  *T = chk_T[ stack_ptr ];
826
827  for(i=0; i<NVAR*S; i++) {
828    Ystage[i] = chk_Y[stack_ptr][i];
829    K[i]      = chk_K[stack_ptr][i];
830  }
831
832  if (SaveLU) {
833#ifdef FULL_ALGEBRA
834    int j;
835    for(i=0; i<NVAR; i++) {
836      for(j=0; j<NVAR; j++)
837        E[(j*NVAR)+i] = chk_J[stack_ptr][j][i];
838      P[i] = chk_P[stack_ptr][i];
839    }
840#else
841    for(i=0; i<LU_NONZERO; i++)
842      E[i] = chk_J[stack_ptr][i];
843#endif
844  }
845
846  stack_ptr--;
847
848} /* End of ros_DPop */
849
850/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
851void ros_CPush( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], 
852                KPP_REAL d2Y[] ) {
853/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
854~~~> Saves the next trajectory snapshot for discrete adjoints
855~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
856   
857  int i;
858 
859  stack_ptr++;
860  if ( stack_ptr > bufsize ) {
861    printf( "Push failed: buffer overflow" );
862    exit(0);
863  }
864  chk_H[ stack_ptr ] = H;
865  chk_T[ stack_ptr ] = T;
866
867  for(i = 0; i< NVAR; i++ ) {
868    chk_Y[stack_ptr][i] = Y[i];
869    chk_dY[stack_ptr][i] = dY[i];
870    chk_d2Y[stack_ptr][i] = d2Y[i];
871  }
872} /* End of ros_CPush */
873
874/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
875void ros_CPop( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL dY[], 
876               KPP_REAL d2Y[] ) {
877/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
878~~~> Retrieves the next trajectory snapshot for discrete adjoints
879~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
880   
881  int i;
882
883  if ( stack_ptr <= 0 ) {
884    printf( "Pop failed: empty buffer" );
885    exit(0);
886  }
887  H = chk_H[ stack_ptr ];
888  T = chk_T[ stack_ptr ];
889
890  for(i=0; i<NVAR; i++) {
891    Y[i] = chk_Y[stack_ptr][i];
892    dY[i] = chk_dY[stack_ptr][i];
893    d2Y[i] = chk_d2Y[stack_ptr][i];
894  }
895  stack_ptr--;
896} /* End of ros_CPop */
897
898/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
899int ros_ErrorMsg( int Code, KPP_REAL T, KPP_REAL H) {
900/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
901    Handles all error messages
902~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
903
904  int IERR = Code;
905  printf( "Forced exit from RosenbrockADJ due to the following error:");
906 
907  switch (Code) {
908    case -1:
909      printf( "--> Improper value for maximal no of steps" );
910      break;
911    case -2:
912      printf( "--> Selected RosenbrockADJ method not implemented" );
913      break;
914    case -3:
915      printf( "--> Hmin/Hmax/Hstart must be positive" );
916      break;
917    case -4:
918      printf( "--> FacMin/FacMax/FacRej must be positive" );
919      break;
920    case -5:
921      printf( "--> Improper tolerance values" );
922      break;
923    case -6:
924      printf( "--> No of steps exceeds maximum buffer bound" );
925      break;
926    case -7:
927      printf( "--> Step size too small: T + 10*H = T or H < Roundoff" );
928      break;
929    case -8:
930      printf( "--> Matrix is repeatedly singular" );
931      break;
932    case -9:
933      printf( "--> Improper type of adjoint selected" );
934      break;
935    default:
936      printf( "Unknown Error code: %d", Code );
937  } /* End of switch */
938
939  return IERR;
940} /* End of ros_ErrorMsg */
941
942/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
943int ros_FwdInt ( KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL T, 
944                 KPP_REAL AbsTol[], KPP_REAL RelTol[], int AdjointType, 
945                 KPP_REAL Hmin, KPP_REAL Hstart, KPP_REAL Hmax, 
946                 KPP_REAL Roundoff, int ISTATUS[], int Max_no_steps, 
947                 KPP_REAL RSTATUS[], int Autonomous, int VectorTol, 
948                 KPP_REAL FacMax, KPP_REAL FacMin, KPP_REAL FacSafe, 
949                 KPP_REAL FacRej, int SaveLU ) {
950/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
951   Template for the implementation of a generic RosenbrockADJ method
952      defined by ros_S (no of stages)
953      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
954~~~> Y - Input: the initial condition at Tstart; Output: the solution at T
955~~~> Tstart, Tend - Input: integration interval
956~~~> T - Output: time at which the solution is returned (T=Tend if success)
957~~~> AbsTol, RelTol - Input: tolerances
958~~~> IERR - Output: Error indicator
959~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
960
961/* ~~~~ Local variables */
962  KPP_REAL Ynew[NVAR], Fcn0[NVAR], Fcn[NVAR];
963  KPP_REAL K[NVAR*ros_S], dFdT[NVAR];
964  KPP_REAL *Ystage = NULL; /* Array pointer */
965#ifdef FULL_ALGEBRA
966  KPP_REAL Jac0[NVAR][NVAR],  Ghimj[NVAR][NVAR];
967#else
968  KPP_REAL Jac0[LU_NONZERO], Ghimj[LU_NONZERO];
969#endif
970  KPP_REAL H, Hnew, HC, HG, Fac, Tau;
971  KPP_REAL Err, Yerr[NVAR];
972  int Pivot[NVAR], Direction, ioffset, i, j=0, istage;
973  int RejectLastH, RejectMoreH, Singular; /* Boolean Values */
974/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
975
976/*~~~>  Allocate stage vector buffer if needed */
977  if (AdjointType == Adj_discrete) {
978    Ystage = (KPP_REAL*) malloc(NVAR*ros_S*sizeof(KPP_REAL));
979    /* Uninitialized Ystage may lead to NaN on some compilers */
980    if (Ystage == NULL) {
981      printf( "Allocation of Ystage failed" );
982      exit(0);
983    }
984  }
985
986/*~~~>  Initial preparations */
987  T = Tstart;
988  RSTATUS[Nhexit] = ZERO;
989  H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) );
990  if (ABS(H) <= ((KPP_REAL)10.0)*Roundoff)
991    H = DeltaMin;
992
993  if (Tend >= Tstart)
994    Direction = 1;
995  else
996    Direction = -1;
997
998  H = Direction*H;
999 
1000  RejectLastH = FALSE;
1001  RejectMoreH = FALSE;
1002
1003/*~~~> Time loop begins below */
1004  while ( ((Direction > 0) && ((T-Tend)+Roundoff <= ZERO)) || 
1005          ((Direction < 0) && ((Tend-T)+Roundoff <= ZERO)) ) { /* TimeLoop */
1006
1007    if ( ISTATUS[Nstp] > Max_no_steps )  /* Too many steps */
1008      return ros_ErrorMsg(-6,T,H);
1009
1010    if ( ((T+((KPP_REAL)0.1)*H) == T) || (H <= Roundoff) ) /* Step size
1011                                                            too small */
1012      return ros_ErrorMsg(-7,T,H);
1013
1014/*~~~>  Limit H if necessary to avoid going beyond Tend */
1015    RSTATUS[Nhexit] = H;
1016    H = MIN(H,ABS(Tend-T));
1017
1018/*~~~>   Compute the function at current time */
1019    FunTemplate(T,Y,Fcn0);
1020    ISTATUS[Nfun] = ISTATUS[Nfun] + 1;
1021
1022/*~~~>  Compute the function derivative with respect to T */
1023    if (!Autonomous)
1024      ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT, ISTATUS );
1025
1026/*~~~>   Compute the Jacobian at current time */
1027    JacTemplate(T,Y,Jac0);
1028    ISTATUS[Njac] = ISTATUS[Njac] + 1;
1029
1030/*~~~>  Repeat step calculation until current step accepted */
1031    do {  /* UntilAccepted */
1032
1033      Singular = ros_PrepareMatrix ( H,Direction,ros_Gamma[0],Jac0,Ghimj,Pivot,
1034                                     ISTATUS );
1035
1036      if (Singular) /* More than 5 consecutive failed decompositions */
1037        return ros_ErrorMsg(-8,T,H);
1038
1039/*~~~>   Compute the stages */
1040      for( istage = 0; istage < ros_S; istage++ ) { /* Stage */
1041
1042        /* Current istage offset. Current istage vector is
1043           K(ioffset+1:ioffset+NVAR) */
1044        ioffset = NVAR*istage;
1045
1046        /*For the 1st istage the function has been computed previously*/
1047        if ( istage == 0 ) {
1048          WCOPY(NVAR,Fcn0,1,Fcn,1);
1049          if (AdjointType == Adj_discrete) { /* Save stage solution */
1050            for(i=0; i<NVAR; i++)
1051              Ystage[i] = Y[i];
1052            WCOPY(NVAR,Y,1,Ynew,1);
1053          }
1054        }
1055        /* istage>0 and a new function evaluation is needed at the
1056           current istage */
1057        else if ( ros_NewF[istage] ) {
1058          WCOPY(NVAR,Y,1,Ynew,1);
1059          for ( j = 0; j < istage; j++ ) {
1060            WAXPY( NVAR,ros_A[(istage)*(istage-1)/2+j], 
1061                   &K[NVAR*j],1,Ynew,1 );
1062          }
1063          Tau = T + ros_Alpha[istage]*Direction*H;
1064          FunTemplate(Tau,Ynew,Fcn);
1065          ISTATUS[Nfun] = ISTATUS[Nfun] + 1;
1066        } /* if istage == 1 elseif ros_NewF[istage] */
1067
1068        /* Save stage solution every time even if ynew is not updated */
1069        if ( ( istage > 0 ) && (AdjointType == Adj_discrete) ) {
1070          for(i=0; i<NVAR; i++)
1071            Ystage[ioffset+i] = Ynew[i];
1072        }
1073        WCOPY(NVAR,Fcn,1,&K[ioffset],1);
1074        for( j = 0; j < istage; j++ ) {
1075          HC = ros_C[(istage)*(istage-1)/2+j]/(Direction*H);
1076          WAXPY(NVAR,HC,&K[NVAR*j],1,&K[ioffset],1);
1077        }
1078        if (( !Autonomous) && (ros_Gamma[istage] != ZERO)) {
1079          HG = Direction*H*ros_Gamma[istage];
1080          WAXPY(NVAR,HG,dFdT,1,&K[ioffset],1);
1081        }
1082        ros_Solve('N', Ghimj, Pivot, &K[ioffset], ISTATUS);
1083      } /* End of Stage loop */
1084
1085/*~~~>  Compute the new solution */
1086      WCOPY(NVAR,Y,1,Ynew,1);
1087      for( j=0; j<ros_S; j++ )
1088        WAXPY(NVAR,ros_M[j],&K[NVAR*j],1,Ynew,1);
1089
1090/*~~~>  Compute the error estimation */ 
1091      WSCAL(NVAR,ZERO,Yerr,1);
1092      for( j=0; j<ros_S; j++ )
1093        WAXPY(NVAR,ros_E[j],&K[NVAR*j],1,Yerr,1);
1094      Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol );
1095
1096/*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax */
1097      Fac  = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,(ONE/ros_ELO))));
1098      Hnew = H*Fac;
1099
1100/*~~~>  Check the error magnitude and adjust step size */
1101      ISTATUS[Nstp] = ISTATUS[Nstp] + 1;
1102      if ( (Err <= ONE) || (H <= Hmin) ) {  /*~~~> Accept step */
1103        ISTATUS[Nacc]++;
1104        if (AdjointType == Adj_discrete) { /* Save current state */
1105          ros_DPush( ros_S, T, H, Ystage, K, Ghimj, Pivot, SaveLU );
1106        }
1107        else if ( (AdjointType == Adj_continuous) || 
1108                  (AdjointType == Adj_simple_continuous) ) {
1109#ifdef FULL_ALGEBRA
1110          K = MATMUL(Jac0,Fcn0);
1111#else
1112          Jac_SP_Vec( Jac0, Fcn0, &K[0] );
1113#endif
1114          if ( !Autonomous)
1115            WAXPY(NVAR,ONE,dFdT,1,&K[0],1);
1116          ros_CPush( T, H, Y, Fcn0, &K[0] );
1117        }
1118        WCOPY(NVAR,Ynew,1,Y,1);
1119        T = T + Direction*H;
1120        Hnew = MAX(Hmin,MIN(Hnew,Hmax));
1121        if (RejectLastH) { /* No step size increase after a
1122                             rejected step */
1123          Hnew = MIN(Hnew,H);
1124        }
1125        RSTATUS[Nhexit] = H;
1126        RSTATUS[Nhnew]  = Hnew;
1127        RSTATUS[Ntexit] = T;
1128        RejectLastH = FALSE;
1129        RejectMoreH = FALSE;
1130        H = Hnew;
1131        break; /* UntilAccepted - EXIT THE LOOP: WHILE STEP NOT ACCEPTED */
1132      }
1133
1134      else { /*~~~> Reject step */
1135        if (RejectMoreH)
1136          Hnew = H*FacRej;
1137        RejectMoreH = RejectLastH;
1138        RejectLastH = TRUE;
1139        H = Hnew;
1140        if (ISTATUS[Nacc] >= 1)
1141          ISTATUS[Nrej]++;
1142      } /* End if else - Err <= 1 */
1143   
1144    } while(1); /* End of UntilAccepted do loop */
1145  } /* End of TimeLoop */
1146   
1147/*~~~> Save last state: only needed for continuous adjoint */
1148  if ( (AdjointType == Adj_continuous) || 
1149       (AdjointType == Adj_simple_continuous) ) {
1150    FunTemplate(T,Y,Fcn0);
1151    ISTATUS[Nfun]++;
1152    JacTemplate(T,Y,Jac0);
1153    ISTATUS[Njac]++;
1154#ifdef FULL_ALGEBRA
1155    K = MATMUL(Jac0,Fcn0);
1156#else
1157    Jac_SP_Vec( Jac0, Fcn0, &K[0] );
1158#endif
1159    if (!Autonomous) {
1160      ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT, ISTATUS );
1161      WAXPY(NVAR,ONE,dFdT,1,&K[0],1);
1162    }
1163    ros_CPush( T, H, Y, Fcn0, &K[0] );
1164/*~~~> Deallocate stage buffer: only needed for discrete adjoint */
1165  }
1166  else if (AdjointType == Adj_discrete) {
1167    free(Ystage);
1168  }
1169
1170/*~~~> Succesful exit */
1171  return 1;  /*~~~> The integration was successful */
1172} /* End of ros_FwdInt */
1173
1174/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1175int ros_DadjInt ( int NADJ, KPP_REAL Lambda[][NVAR], KPP_REAL Tstart, 
1176                  KPP_REAL Tend, KPP_REAL T, int SaveLU, int ISTATUS[], 
1177                  KPP_REAL Roundoff, int Autonomous) {
1178/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1179   Template for the implementation of a generic RosenbrockSOA method
1180      defined by ros_S (no of stages)
1181      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1182!~~~> NADJ - Input: the initial condition at Tstart; Output: the solution at T
1183!~~~> Lambda[NADJ][NVAR] -  First order adjoint
1184!~~~> Tstart, Tend - Input: integration interval
1185!~~~> T - Output: time at which the solution is returned (T=Tend if success)
1186~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1187
1188/*~~~~ Local variables */
1189  KPP_REAL Ystage[NVAR*ros_S], K[NVAR*ros_S];
1190  KPP_REAL U[NADJ][NVAR*ros_S], V[NADJ][NVAR*ros_S];
1191#ifdef FULL_ALGEBRA
1192  KPP_REAL Jac[NVAR][NVAR], dJdT[NVAR][NVAR], Ghimj[NVAR][NVAR];
1193#else
1194  KPP_REAL Jac[LU_NONZERO], dJdT[LU_NONZERO], Ghimj[LU_NONZERO];
1195#endif
1196  KPP_REAL Hes0[NHESS];
1197  KPP_REAL Tmp[NVAR], Tmp2[NVAR];
1198  KPP_REAL H=0.0, HC, HA, Tau;
1199  int Pivot[NVAR], Direction;
1200  int i, j, m, istage, istart, jstart;
1201/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1202
1203  if (Tend  >=  Tstart)
1204    Direction = 1;
1205  else
1206    Direction = -1;
1207 
1208  /*~~~> Time loop begins below */
1209  while ( stack_ptr > 0 ) { /* TimeLoop */
1210
1211    /*~~~>  Recover checkpoints for stage values and vectors */
1212    ros_DPop( ros_S, &T, &H, &Ystage[0], &K[0], &Ghimj[0], &Pivot[0], SaveLU );
1213
1214    /*~~~>    Compute LU decomposition */
1215    if (!SaveLU) {
1216      JacTemplate(T,&Ystage[0],Ghimj);
1217      ISTATUS[Njac] = ISTATUS[Njac] + 1;
1218      Tau = ONE/(Direction*H*ros_Gamma[0]);
1219#ifdef FULL_ALGEBRA
1220      for(j=0; j<NVAR; j++) {
1221        for(i=0; i<NVAR; i++)
1222          Ghimj[i][j] = -Ghimj[i][j];
1223      }
1224      for(i=0; i<NVAR; i++)
1225        Ghimj[i][i] = Ghimj[i][i]+Tau;
1226#else
1227      WSCAL(LU_NONZERO,(-ONE),Ghimj,1);
1228      for (i=0; i<NVAR; i++)
1229        Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+Tau;
1230#endif
1231      ros_Decomp(Ghimj, Pivot, &j, ISTATUS);
1232    }
1233
1234/*~~~>   Compute Hessian at the beginning of the interval */
1235    HessTemplate(T,&Ystage[0],Hes0);
1236
1237/*~~~>   Compute the stages */
1238    for (istage = ros_S - 1; istage >= 0; istage--) { /* Stage loop */
1239
1240      /*~~~> Current istage first entry */
1241      istart = NVAR*istage;
1242         
1243      /*~~~> Compute U */
1244      for (m = 0; m<NADJ; m++) {
1245        WCOPY(NVAR,&Lambda[m][0],1,&U[m][istart],1);
1246        WSCAL(NVAR,ros_M[istage],&U[m][istart],1);
1247      } /* m=0:NADJ-1 */
1248      for (j = istage+1; j < ros_S; j++) {
1249        jstart = NVAR*j;
1250        HA = ros_A[j*(j-1)/2+istage];
1251        HC = ros_C[j*(j-1)/2+istage]/(Direction*H);
1252        for ( m = 0; m < NADJ; m++ ) {
1253          WAXPY(NVAR,HA,&V[m][jstart],1,&U[m][istart],1);
1254          WAXPY(NVAR,HC,&U[m][jstart],1,&U[m][istart],1);
1255        } /* m=0:NADJ-1 */
1256      }
1257      for ( m = 0; m < NADJ; m++ )
1258        ros_Solve('T', Ghimj, Pivot, &U[m][istart], ISTATUS); /* m=1:NADJ-1 */
1259
1260      /*~~~> Compute V */
1261      Tau = T + ros_Alpha[istage]*Direction*H;
1262      JacTemplate(Tau,&Ystage[istart],Jac);
1263      ISTATUS[Njac]++;
1264      for ( m = 0; m < NADJ; m++ ) {
1265#ifdef FULL_ALGEBRA
1266        for (i=istart; i < istart+NVAR-1; i++ )
1267          V[[m][i] = MATMUL(TRANSPOSE(Jac),U[m][i]];
1268#else
1269        JacTR_SP_Vec(Jac,&U[m][istart],&V[m][istart]); 
1270#endif
1271      } /* m=0:NADJ-1 */
1272    } /*End of Stage loop */
1273
1274    if (!Autonomous)
1275/*~~~>  Compute the Jacobian derivative with respect to T.
1276        Last "Jac" computed for stage 1 */
1277      ros_JacTimeDerivative ( T, Roundoff, &Ystage[0], Jac, dJdT, ISTATUS );
1278
1279/*~~~>  Compute the new solution */
1280    /*~~~>  Compute Lambda */ 
1281    for( istage = 0; istage < ros_S; istage++ ) {
1282      istart = NVAR*istage;
1283      for (m = 0; m < NADJ; m++) {
1284        /* Add V_i */
1285        WAXPY(NVAR,ONE,&V[m][istart],1,&Lambda[m][0],1);
1286        /* Add (H0xK_i)^T * U_i */
1287        HessTR_Vec ( Hes0, &U[m][istart], &K[istart], Tmp );
1288        WAXPY(NVAR,ONE,Tmp,1,&Lambda[m][0],1);
1289      } /* m=0:NADJ-1 */
1290    }
1291
1292    /* Add H * dJac_dT_0^T * \sum(gamma_i U_i) */
1293    /* Tmp holds sum gamma_i U_i */
1294    if (!Autonomous) {
1295      for( m = 0; m < NADJ; m++ ) {
1296        for(i=0; i<NVAR; i++)
1297          Tmp[i] = ZERO;
1298        for( istage = 0; istage < ros_S; istage++ ) {
1299          istart = NVAR*istage;
1300          WAXPY(NVAR,ros_Gamma[istage],&U[m][istart],1,Tmp,1);
1301        }
1302#ifdef FULL_ALGEBRA
1303        Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp);
1304#else
1305        JacTR_SP_Vec(dJdT,Tmp,Tmp2);
1306#endif
1307        WAXPY(NVAR,H,Tmp2,1,&Lambda[m][0],1);
1308      } /* m=0:NADJ-1 */
1309    } /* .NOT.Autonomous */
1310  } /* End of TimeLoop */
1311
1312  /*~~~> Save last state */
1313  /*~~~> Succesful exit */
1314  return 1;  /*~~~> The integration was successful */
1315
1316/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1317} /* End of ros_DadjInt */
1318/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1319
1320/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1321int ros_CadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart, KPP_REAL Tend,
1322                  KPP_REAL T, KPP_REAL AbsTol_adj[][NVAR], 
1323                  KPP_REAL RelTol_adj[][NVAR], KPP_REAL RSTATUS[], 
1324                  KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, 
1325                  KPP_REAL Roundoff, int Max_no_steps, int Autonomous, 
1326                  int VectorTol, KPP_REAL FacMax, KPP_REAL FacMin, 
1327                  KPP_REAL FacSafe, KPP_REAL FacRej, int ISTATUS[] ) {
1328/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1329   Template for the implementation of a generic RosenbrockADJ method
1330      defined by ros_S (no of stages)
1331      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1332~~~> NADJ, Y[NADJ][NVAR] - Input: the initial condition at Tstart;
1333                           Output: the solution at T
1334~~~> Tstart, Tend - Input: integration interval
1335~~~> AbsTol_adj[NADJ][NVAR], RelTol_adj[NADJ][NVAR] - Input: adjoint tolerances
1336~~~> T - Output: time at which the solution is returned (T=Tend if success)
1337~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1338
1339/*~~~~ Local variables */
1340  KPP_REAL Y0[NVAR];
1341  KPP_REAL Ynew[NADJ][NVAR], Fcn0[NADJ][NVAR], Fcn[NADJ][NVAR];
1342  KPP_REAL K[NADJ][NVAR*ros_S], dFdT[NADJ][NVAR];
1343#ifdef FULL_ALGEBRA
1344  KPP_REAL Jac0[NVAR][NVAR], Ghimj[NVAR][NVAR], Jac[NVAR][NVAR], 
1345         dJdT[NVAR][NVAR];
1346#else
1347  KPP_REAL Jac0[LU_NONZERO], Ghimj[LU_NONZERO], Jac[LU_NONZERO], 
1348         dJdT[LU_NONZERO];
1349#endif
1350  KPP_REAL H, Hnew, HC, HG, Fac, Tau; 
1351  KPP_REAL Err, Yerr[NADJ][NVAR];
1352  int Pivot[NVAR], Direction, ioffset, j, istage, iadj;
1353  int RejectLastH, RejectMoreH, Singular; /* Boolean values */
1354/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1355
1356/*~~~>  Initial preparations */
1357  T = Tstart;
1358  RSTATUS[Nhexit] = (KPP_REAL)0.0;
1359  H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) );
1360  if (ABS(H) <= ((KPP_REAL)10.0)*Roundoff) 
1361    H = DeltaMin;
1362 
1363  if (Tend  >=  Tstart)
1364    Direction = 1;
1365  else
1366    Direction = -1;
1367 
1368  H = Direction*H;
1369  RejectLastH = FALSE;
1370  RejectMoreH = FALSE;
1371
1372/*~~~> Time loop begins below */
1373  while ( ((Direction > 0) && ((T-Tend)+Roundoff <= ZERO))
1374          || ((Direction < 0) && ((Tend-T)+Roundoff <= ZERO)) ) { /*TimeLoop*/
1375
1376    if ( ISTATUS[Nstp] > Max_no_steps )  /* Too many steps */
1377      return ros_ErrorMsg(-6,T,H);
1378
1379    /* Step size too small */
1380    if ( ((T+((KPP_REAL)0.1)*H) == T) || (H <= Roundoff) )
1381      return ros_ErrorMsg(-7,T,H);
1382
1383/*~~~>  Limit H if necessary to avoid going beyond Tend */
1384    RSTATUS[Nhexit] = H;
1385    H = MIN(H,ABS(Tend-T));
1386
1387/*~~~>   Interpolate forward solution */
1388    ros_cadj_Y( T, Y0 );
1389/*~~~>   Compute the Jacobian at current time */
1390    JacTemplate(T, Y0, Jac0);
1391    ISTATUS[Njac]++;
1392
1393/*~~~>  Compute the function derivative with respect to T */
1394    if (!Autonomous) {
1395      ros_JacTimeDerivative ( T, Roundoff, Y0, Jac0, dJdT, ISTATUS );
1396      for (iadj = 0; iadj < NADJ; iadj++) {
1397#ifdef FULL_ALGEBRA
1398        for (i=0; i<NVAR; i++)
1399          dFdT[iadj][i] = MATMUL(TRANSPOSE(dJdT),Y[iadj][i]);
1400#else
1401        JacTR_SP_Vec(dJdT,&Y[iadj][0],&dFdT[iadj][0]);
1402#endif
1403        WSCAL(NVAR,(-ONE),&dFdT[iadj][0],1);
1404      } /* End for loop */
1405    } /* End if */
1406
1407/*~~~>  Ydot = -J^T*Y */
1408#ifdef FULL_ALGEBRA
1409    int i;
1410    for(i=0; i<NVAR; i++) {
1411      for(j=0; j<NVAR; j++)
1412        Jac0[i][j] = -Jac0[i][j];
1413    }
1414#else
1415    WSCAL(LU_NONZERO,(-ONE),Jac0,1);
1416#endif
1417
1418    for ( iadj = 0; iadj < NADJ; iadj++ ) {
1419#ifdef FULL_ALGEBRA
1420      int i;
1421      for(i=0; i<NVAR; i++)
1422        Fcn0[iadj][i] = MATMUL(TRANSPOSE(Jac0),Y[iadj][i]);
1423#else
1424      JacTR_SP_Vec(Jac0,&Y[iadj][0],&Fcn0[iadj][0]);
1425#endif
1426    }
1427
1428/*~~~>  Repeat step calculation until current step accepted */
1429    do { /* UntilAccepted */
1430
1431      Singular = ros_PrepareMatrix(H,Direction,ros_Gamma[0], Jac0,Ghimj,Pivot,
1432                                   ISTATUS);
1433      if (Singular) /* More than 5 consecutive failed decompositions */
1434        return ros_ErrorMsg(-8,T,H);
1435
1436/*~~~>   Compute the stages */
1437      for ( istage = 0; istage < ros_S; istage++ ) { /* Stage loop */
1438
1439        /* Current istage offset. Current istage vector
1440           is K[ioffset+1:ioffset+NVAR] */
1441        ioffset = NVAR*(istage-1);
1442
1443        /* For the 1st istage the function has been computed previously */
1444        if ( istage == 0 ) {
1445          for ( iadj = 0; iadj < NADJ; iadj++ )
1446            WCOPY(NVAR,&Fcn0[iadj][0],1,&Fcn[iadj][0],1);
1447
1448          /* istage>0 and a new function evaluation is needed at
1449             the current istage */
1450        }
1451        else if ( ros_NewF[istage] ) {
1452          WCOPY(NVAR*NADJ,&Y[0][0],1,&Ynew[0][0],1);
1453          for (j = 0; j < istage-1; j++) {
1454            for ( iadj = 0; iadj < NADJ; iadj++ )
1455              WAXPY( NVAR,ros_A[(istage-1)*(istage-2)/2+j],
1456                     &K[iadj][NVAR*(j-1)+1],1,&Ynew[iadj][0],1);
1457          } /* End for loop */
1458          Tau = T + ros_Alpha[istage]*Direction*H;
1459          ros_cadj_Y( Tau, Y0 );
1460          JacTemplate(Tau, Y0, Jac);
1461          ISTATUS[Njac]++;
1462
1463#ifdef FULL_ALGEBRA
1464          for(i=0; i<NVAR; i++) {
1465            for(j=0; j<NVAR; j++)
1466              Jac[i][j] = -Jac[i][j];
1467          }
1468#else
1469          WSCAL(LU_NONZERO,(-ONE),Jac,1);
1470#endif
1471
1472          for ( iadj = 0; iadj < NADJ; iadj++ ) {
1473#ifdef FULL_ALGEBRA
1474            for(i=0; i<NVAR; i++)
1475              Fcn[iadj][i] = MATMUL(TRANSPOSE(Jac),Ynew[iadj][i]);
1476#else
1477            JacTR_SP_Vec(Jac,&Ynew[iadj][0],&Fcn[iadj][0]);
1478#endif
1479          } /* End for loop */
1480        } /*  if istage == 1 elseif ros_NewF(istage) */
1481
1482        for ( iadj = 0; iadj < NADJ; iadj++ )
1483          WCOPY(NVAR,&Fcn[iadj][0],1,&K[iadj][ioffset+1],1);
1484        for ( j = 0; j < istage-1; j++ ) {
1485          HC = ros_C[(istage-1)*(istage-2)/2+j]/(Direction*H);
1486          for ( iadj = 0; iadj < NADJ; iadj++ )
1487            WAXPY(NVAR,HC,&K[iadj][NVAR*(j-1)+1],1,&K[iadj][ioffset+1],1);
1488        } /* End for loop */
1489        if ((!Autonomous) && (ros_Gamma[istage] != ZERO)) {
1490          HG = Direction*H*ros_Gamma[istage];
1491          for ( iadj = 0; iadj < NADJ; iadj++ )
1492            WAXPY(NVAR,HG,&dFdT[iadj][0],1,&K[iadj][ioffset+1],1);
1493        } /* End if */
1494        for ( iadj = 0; iadj < NADJ; iadj++ )
1495          ros_Solve('T', Ghimj, Pivot, &K[iadj][ioffset+1], ISTATUS);
1496
1497      } /* End of Stage loop */
1498
1499/*~~~>  Compute the new solution */
1500      for ( iadj = 0; iadj < NADJ; iadj++ ) {
1501        WCOPY(NVAR,&Y[iadj][0],1,&Ynew[iadj][0],1);
1502        for ( j=0; j<ros_S; j++ )
1503          WAXPY(NVAR,ros_M[j],&K[iadj][NVAR*(j-1)+1],1,&Ynew[iadj][0],1);
1504      } /* End for loop */
1505
1506/*~~~>  Compute the error estimation */
1507      WSCAL(NVAR*NADJ,ZERO,&Yerr[0][0],1);
1508      for ( j=0; j<ros_S; j++ ) {
1509        for ( iadj = 0; iadj < NADJ; iadj++ )
1510          WAXPY(NVAR,ros_E[j],&K[iadj][NVAR*(j-1)+1],1,&Yerr[iadj][0],1);
1511      } /* End for loop */
1512
1513/*~~~> Max error among all adjoint components */
1514      iadj = 1;
1515      Err = ros_ErrorNorm ( &Y[iadj][0], &Ynew[iadj][0], &Yerr[iadj][0],
1516                            &AbsTol_adj[iadj][0], &RelTol_adj[iadj][0], 
1517                            VectorTol );
1518
1519/*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax */
1520      Fac  = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,(ONE/ros_ELO))));
1521      Hnew = H*Fac;
1522
1523/*~~~>  Check the error magnitude and adjust step size */
1524      /* ISTATUS[Nstp] = ISTATUS[Nstp] + 1 */
1525      if ( (Err <= ONE) || (H <= Hmin) ) {  /*~~~> Accept step */
1526        ISTATUS[Nacc] = ISTATUS[Nacc] + 1;
1527        WCOPY(NVAR*NADJ,&Ynew[0][0],1,&Y[0][0],1);
1528        T = T + Direction*H;
1529        Hnew = MAX(Hmin,MIN(Hnew,Hmax));
1530        if (RejectLastH) /* No step size increase after a rejected step */
1531          Hnew = MIN(Hnew,H);
1532        RSTATUS[Nhexit] = H;
1533        RSTATUS[Nhnew] = Hnew;
1534        RSTATUS[Ntexit] = T;
1535        RejectLastH = FALSE;
1536        RejectMoreH = FALSE;
1537        H = Hnew;
1538        break; /* UntilAccepted - EXIT THE LOOP: WHILE STEP NOT ACCEPTED */
1539      }
1540      else {         /*~~~> Reject step */
1541        if (RejectMoreH)
1542          Hnew = H*FacRej;
1543        RejectMoreH = RejectLastH;
1544        RejectLastH = TRUE;
1545        H = Hnew;
1546        if (ISTATUS[Nacc] >= 1)
1547          ISTATUS[Nrej]++;
1548      } /* Err <= 1 */
1549     
1550    } while(1); /* End of UntilAccepted do loop */
1551
1552  } /* End of TimeLoop */
1553
1554  /*~~~> Succesful exit */
1555  return 1;  /*~~~> The integration was successful */
1556} /* End of ros_CadjInt */
1557
1558/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1559int ros_SimpleCadjInt ( int NADJ, KPP_REAL Y[][NVAR], KPP_REAL Tstart,
1560                         KPP_REAL Tend, KPP_REAL T, int ISTATUS[], 
1561                        int Autonomous, KPP_REAL Roundoff ) {
1562/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1563   Template for the implementation of a generic RosenbrockADJ method
1564      defined by ros_S (no of stages)
1565      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1566~~~> NADJ, Y[NADJ][NVAR] - Input: the initial condition at Tstart;
1567                           Output: the solution at T
1568~~~> Tstart, Tend - Input: integration interval
1569~~~> T - Output: time at which the solution is returned (T=Tend if success)
1570~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1571
1572/*~~~~ Local variables */
1573   KPP_REAL Y0[NVAR];
1574   KPP_REAL Ynew[NADJ][NVAR], Fcn0[NADJ][NVAR], Fcn[NADJ][NVAR];
1575   KPP_REAL K[NADJ][NVAR*ros_S], dFdT[NADJ][NVAR];
1576#ifdef FULL_ALGEBRA
1577   KPP_REAL Jac0[NVAR][NVAR], Ghimj[NVAR][NVAR], Jac[NVAR][NVAR], 
1578          dJdT[NVAR][NVAR];
1579#else
1580   KPP_REAL Jac0[LU_NONZERO], Ghimj[LU_NONZERO], Jac[LU_NONZERO], 
1581          dJdT[LU_NONZERO];
1582#endif
1583   KPP_REAL H, HC, HG, Tau; 
1584   KPP_REAL ghinv;
1585   int Pivot[NVAR], Direction, ioffset, i, j, istage, iadj;
1586   int istack;
1587
1588/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1589
1590/*~~~>  INITIAL PREPARATIONS */
1591  if (Tend  >=  Tstart)
1592    Direction = -1;
1593  else
1594    Direction = 1;
1595 
1596/*~~~> Time loop begins below */
1597  for( istack = stack_ptr; istack >= 1; istack-- ) { /* TimeLoop */
1598    T = chk_T[istack];
1599    H = chk_H[istack-1];
1600    for(i=0; i<NVAR; i++)
1601      Y0[i] = chk_Y[istack][i];
1602
1603/*~~~>   Compute the Jacobian at current time */
1604    JacTemplate(T, Y0, Jac0);
1605    ISTATUS[Njac] = ISTATUS[Njac] + 1;
1606
1607/*~~~>  Compute the function derivative with respect to T */
1608    if (!Autonomous) {
1609      ros_JacTimeDerivative ( T, Roundoff, Y0, Jac0, dJdT, ISTATUS );
1610      for ( iadj = 0; iadj < NADJ; iadj++ ) {
1611#ifdef FULL_ALGEBRA
1612        for( i=0; i<NVAR; i++)
1613          dFdT[iadj][i] = MATMUL(TRANSPOSE(dJdT),Y[iadj][i]);
1614#else
1615        JacTR_SP_Vec(dJdT,&Y[iadj][0],&dFdT[iadj][0]);
1616#endif
1617        WSCAL(NVAR,(-ONE),&dFdT[iadj][0],1);
1618      }
1619    }
1620
1621/*~~~>  Ydot = -J^T*Y */
1622#ifdef FULL_ALGEBRA
1623    for(i=0; i<NVAR; i++) {
1624      for(j=0; j<NVAR; j++)
1625        Jac0[i][j] = -Jac0[i][j];
1626    }
1627#else
1628    WSCAL(LU_NONZERO,(-ONE),Jac0,1);
1629#endif
1630   
1631    for(iadj=0; iadj<NADJ; iadj++) {
1632#ifdef FULL_ALGEBRA
1633      for(i=0; i<NVAR; i++)
1634        Fcn0[iadj][i] = MATMUL(TRANSPOSE(Jac0),Y[iadj][i]);
1635#else
1636      JacTR_SP_Vec(Jac0,&Y[iadj][0],&Fcn0[iadj][0]);
1637#endif
1638    }
1639   
1640/*~~~>    Construct Ghimj = 1/(H*ham) - Jac0 */
1641    ghinv = ONE/(Direction*H*ros_Gamma[0]);
1642#ifdef FULL_ALGEBRA
1643    for(i=0; i<NVAR; i++) {
1644      for(j=0; j<NVAR; j++)
1645        Ghimj[i][j] = -Jac0[i][j];
1646    }
1647    for(i=0; i<NVAR; i++)
1648      Ghimj[i][i] = Ghimj[i][i]+ghinv;
1649#else
1650    WCOPY(LU_NONZERO,Jac0,1,Ghimj,1);
1651    WSCAL(LU_NONZERO,(-ONE),Ghimj,1);
1652    for(i=0; i<NVAR; i++)
1653      Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+ghinv;
1654#endif
1655
1656/*~~~>    Compute LU decomposition */
1657    ros_Decomp( Ghimj, Pivot, &j, ISTATUS );
1658    if (j != 0) {
1659      ros_ErrorMsg(-8,T,H);
1660      printf( "The matrix is singular !");
1661      exit(0);
1662    }
1663
1664/*~~~>   Compute the stages */
1665    for(istage=0; istage<ros_S; istage++) { /* Stage */
1666      /* Current istage offset. Current istage vector
1667         is K(ioffset+1:ioffset+NVAR) */
1668      ioffset = NVAR*istage;
1669
1670      /* For the 1st istage the function has been computed previously */
1671      if ( istage == 0 ) {
1672        for(iadj=0; iadj<NADJ; iadj++)
1673          WCOPY(NVAR,&Fcn0[iadj][0],1,&Fcn[iadj][0],1);
1674      }
1675      /* istage>=1 and a new function evaluation is needed
1676         at the current istage */
1677      else if ( ros_NewF[istage] ) {
1678        WCOPY(NVAR*NADJ,&Y[0][0],1,&Ynew[0][0],1);
1679        for(j=0; j<istage; j++) {
1680          for(iadj=0; iadj<NADJ; iadj++)
1681            WAXPY(NVAR,ros_A[istage*(istage-1)/2+j], &K[iadj][NVAR*j],1,
1682                  &Ynew[iadj][0],1);
1683        }
1684       
1685        Tau = T + ros_Alpha[istage]*Direction*H;
1686        for(i=0; i<NVAR; i++)
1687          ros_Hermite3( chk_T[istack-1], chk_T[istack], Tau,
1688                        &chk_Y[istack-1][i], &chk_Y[istack][i],
1689                        &chk_dY[istack-1][i], &chk_dY[istack][i], Y0 );
1690        JacTemplate(Tau, Y0, Jac);
1691        ISTATUS[Njac]++;
1692
1693#ifdef FULL_ALGEBRA
1694        for(i=0; i<NVAR; i++) {
1695          for(j=0; j<NVAR; j++)
1696            Jac[i][j] = -Jac[i][j];
1697        }
1698#else
1699        WSCAL(LU_NONZERO,(-ONE),Jac,1);
1700#endif
1701
1702        for(iadj=0; iadj<NADJ; iadj++) {
1703#ifdef FULL_ALGEBRA
1704          for(i=0; i<NVAR; i++)
1705            Fcn[iadj][i] = MATMUL(TRANSPOSE(Jac),Ynew[iadj][i]);
1706#else
1707          JacTR_SP_Vec(Jac,&Ynew[iadj][0],&Fcn[iadj][0]);
1708#endif
1709        }
1710      } /* if istage == 1 elseif ros_NewF(istage) */
1711
1712      for(iadj=0; iadj<NADJ; iadj++)
1713        WCOPY(NVAR,&Fcn[iadj][0],1,&K[iadj][ioffset],1);
1714      for(j=0; j<istage-1; j++) {
1715        HC = ros_C[istage*(istage-1)/2+j]/(Direction*H);
1716        for(iadj=0; iadj<NADJ; iadj++)
1717          WAXPY(NVAR,HC,&K[iadj][NVAR*j],1,&K[iadj][ioffset],1);
1718      }
1719      if((!Autonomous) && (ros_Gamma[istage] != ZERO)) {
1720        HG = Direction*H*ros_Gamma[istage];
1721        for(iadj=0; iadj<NADJ; iadj++)
1722          WAXPY(NVAR,HG,&dFdT[iadj][0],1,&K[iadj][ioffset],1);
1723      }
1724      for(iadj=0; iadj<NADJ; iadj++)
1725        ros_Solve('T', Ghimj, Pivot, &K[iadj][ioffset], ISTATUS);
1726    } /* End of Stage loop */
1727   
1728/*~~~>  Compute the new solution */
1729    for(iadj=0; iadj<NADJ; iadj++) {
1730      for(j=0; j<ros_S; j++)
1731        WAXPY(NVAR,ros_M[j],&K[iadj][NVAR*j],1,&Y[iadj][0],1);
1732    }
1733  } /* End of TimeLoop */
1734
1735/*~~~> Succesful exit */
1736  return 1;  /*~~~> The integration was successful */
1737
1738} /* End of ros_SimpleCadjInt */
1739
1740/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1741KPP_REAL ros_ErrorNorm ( KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[],
1742                       KPP_REAL AbsTol[], KPP_REAL RelTol[], int VectorTol ) {
1743/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1744~~~> Computes the "scaled norm" of the error vector Yerr
1745~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1746
1747/* Local variables */
1748  KPP_REAL Err, Scale, Ymax;
1749  int i;
1750
1751  Err = ZERO;
1752  for(i=0; i<NVAR; i++) {
1753    Ymax = MAX(ABS(Y[i]),ABS(Ynew[i]));
1754    if (VectorTol)
1755      Scale = AbsTol[i]+RelTol[i]*Ymax;
1756    else
1757      Scale = AbsTol[0]+RelTol[0]*Ymax;
1758
1759    Err = Err+pow((Yerr[i]/Scale),2);
1760  }
1761  Err  = SQRT(Err/NVAR);
1762
1763  return MAX(Err,(KPP_REAL)1.0e-10);
1764} /* End of ros_ErrorNorm */
1765
1766/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1767void ros_FunTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[],
1768                             KPP_REAL Fcn0[], KPP_REAL dFdT[], int ISTATUS[]) {
1769/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1770~~~> The time partial derivative of the function by finite differences
1771~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1772
1773/*~~~> Local variables */
1774  KPP_REAL Delta;
1775 
1776  Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T));
1777  FunTemplate(T+Delta,Y,dFdT);
1778  ISTATUS[Nfun]++;
1779  WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1);
1780  WSCAL(NVAR,(ONE/Delta),dFdT,1);
1781
1782} /* End of ros_FunTimeDerivative */
1783
1784/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1785void ros_JacTimeDerivative ( KPP_REAL T, KPP_REAL Roundoff, KPP_REAL Y[],
1786                             KPP_REAL Jac0[], KPP_REAL dJdT[], int ISTATUS[]) {
1787/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1788~~~> The time partial derivative of the Jacobian by finite differences
1789~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1790
1791/*~~~> Local variables */
1792  KPP_REAL Delta;
1793
1794  Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T));
1795  JacTemplate(T+Delta,Y,dJdT);
1796  ISTATUS[Njac]++;
1797#ifdef FULL_ALGEBRA
1798  WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1);
1799  WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1);
1800#else
1801  WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1);
1802  WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1);
1803#endif
1804} /* End of ros_JacTimeDerivative */
1805
1806/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1807int ros_PrepareMatrix ( KPP_REAL H, int Direction, KPP_REAL gam,
1808                        KPP_REAL Jac0[], KPP_REAL Ghimj[], int Pivot[], 
1809                        int ISTATUS[] ) {
1810/* --- --- --- --- --- --- --- --- --- --- --- --- ---
1811  Prepares the LHS matrix for stage calculations
1812  1.  Construct Ghimj = 1/(H*gam) - Jac0
1813      "(Gamma H) Inverse Minus Jacobian"
1814  2.  Repeat LU decomposition of Ghimj until successful.
1815       -half the step size if LU decomposition fails and retry
1816       -exit after 5 consecutive fails
1817 --- --- --- --- --- --- --- --- --- --- --- --- --- */
1818
1819/*~~~> Local variables */
1820  int i, ising, Nconsecutive;
1821  int Singular; /* Boolean value */
1822  KPP_REAL ghinv;
1823
1824  Nconsecutive = 0;
1825  Singular = TRUE;
1826
1827  while (Singular) {
1828
1829/*~~~>    Construct Ghimj = 1/(H*gam) - Jac0 */
1830#ifdef FULL_ALGEBRA
1831    WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1);
1832    WSCAL(NVAR*NVAR,(-ONE),Ghimj,1);
1833    ghinv = ONE/(Direction*H*gam);
1834    for(i=0; i<NVAR; i++)
1835      Ghimj[i][i] = Ghimj[i][i]+ghinv;
1836#else
1837    WCOPY(LU_NONZERO,Jac0,1,Ghimj,1);
1838    WSCAL(LU_NONZERO,(-ONE),Ghimj,1);
1839    ghinv = ONE/(Direction*H*gam);
1840    for(i=0; i<NVAR; i++)
1841      Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+ghinv;
1842#endif
1843
1844/*~~~>    Compute LU decomposition */
1845    ros_Decomp( Ghimj, Pivot, &ising, ISTATUS );
1846    if (ising == 0)
1847/*~~~>    If successful done */
1848      Singular = FALSE;
1849   
1850    else { /* ising != 0 */
1851/*~~~>    If unsuccessful half the step size;
1852          if 5 consecutive fails then return */
1853      ISTATUS[Nsng]++;
1854      Nconsecutive++;
1855      Singular = TRUE;
1856      printf( "Warning: LU Decomposition returned ising = %d", ising );
1857      if (Nconsecutive <= 5) /*Less than 5 consecutive failed decompositions*/
1858        H = H*HALF;
1859      else  /* More than 5 consecutive failed decompositions */
1860        return Singular;
1861
1862    }  /* End of ising */
1863  }  /* while Singular */
1864
1865  return Singular;
1866} /* End of ros_PrepareMatrix */
1867
1868/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1869void ros_Decomp( KPP_REAL A[], int Pivot[], int* ising, int ISTATUS[] ) {
1870/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1871  Template for the LU decomposition
1872~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1873
1874#ifdef FULL_ALGEBRA
1875  DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising );
1876#else
1877  *ising = KppDecomp ( A );
1878  Pivot[0] = 1;
1879#endif
1880  ISTATUS[Ndec]++;
1881
1882} /* End of ros_Decomp */
1883
1884/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1885void ros_Solve( char How, KPP_REAL A[], int Pivot[], KPP_REAL b[], 
1886                int ISTATUS[] ) {
1887/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1888  Template for the forward/backward substitution
1889  (using pre-computed LU decomposition)
1890~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1891  switch (How) {
1892    case 'N':
1893#ifdef FULL_ALGEBRA   
1894      DGETRS( 'N', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 );
1895#else
1896      KppSolve( A, b );
1897#endif
1898      break;
1899    case 'T':
1900#ifdef FULL_ALGEBRA
1901      DGETRS( 'T', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 );
1902#else
1903      KppSolveTR( A, b, b );
1904#endif
1905      break;
1906    default:
1907      printf( "Error: unknown argument in ros_Solve: How=%d", How );
1908      exit(0);
1909  }
1910  ISTATUS[Nsol]++;
1911
1912} /* End of ros_Solve */
1913
1914/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1915void ros_cadj_Y( KPP_REAL T, KPP_REAL Y[] ) {
1916/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1917  Finds the solution Y at T by interpolating the stored forward trajectory
1918~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1919
1920/*~~~> Local variables */
1921  int i,j;
1922
1923  if( (T < chk_T[0]) || (T> chk_T[stack_ptr]) ) {
1924    printf( "Cannot locate solution at T = %f", T );
1925    printf( "Stored trajectory is between Tstart = %f", chk_T[0] );
1926    printf( "    and Tend = %f", chk_T[stack_ptr] );
1927    exit(0);
1928  }
1929  for(i=0; i<stack_ptr-1; i++) {
1930    if( (T >= chk_T[i]) &&(T <= chk_T[i+1]) )
1931      exit(0);
1932  }
1933
1934  for(j=0; j<NVAR; j++ )
1935    ros_Hermite3( chk_T[i], chk_T[i+1], T, &chk_Y[i][j], &chk_Y[i+1][j],
1936                  &chk_dY[i][j],  &chk_dY[i+1][j], Y );
1937
1938} /* End of ros_cadj_Y */
1939
1940/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1941void ros_Hermite3( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[],
1942                   KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Y[]) {
1943/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1944  Template for Hermite interpolation of order 5 on the interval [a,b]
1945   P = c(1) + c(2)*(x-a) + ... + c(4)*(x-a)^3
1946   P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb]
1947~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1948
1949/*~~~> Local variables */
1950  KPP_REAL Tau, amb[3], C[4][NVAR];
1951  int i, j;
1952
1953  amb[0] = ((KPP_REAL)1.0)/(a-b);
1954  for(i=1; i<3; i++)
1955    amb[i] = amb[i-1]*amb[0];
1956
1957/* c(1) = ya; */
1958  WCOPY(NVAR,Ya,1,&C[0][0],1);
1959/* c(2) = ja; */
1960  WCOPY(NVAR,Ja,1,&C[1][0],1);
1961/* c(3) = 2/(a-b)*ja + 1/(a-b)*jb - 3/(a - b)^2*ya + 3/(a - b)^2*yb; */
1962  WCOPY(NVAR,Ya,1,&C[2][0],1);
1963  WSCAL(NVAR,-3.0*amb[1],&C[2][0],1);
1964  WAXPY(NVAR,3.0*amb[1],Yb,1,&C[2][0],1);
1965  WAXPY(NVAR,2.0*amb[0],Ja,1,&C[2][0],1);
1966  WAXPY(NVAR,amb[0],Jb,1,&C[2][0],1);
1967/* c(4) =  1/(a-b)^2*ja + 1/(a-b)^2*jb - 2/(a-b)^3*ya + 2/(a-b)^3*yb */
1968  WCOPY(NVAR,Ya,1,&C[3][0],1);
1969  WSCAL(NVAR,-2.0*amb[2],&C[3][0],1);
1970  WAXPY(NVAR,2.0*amb[2],Yb,1,&C[3][0],1);
1971  WAXPY(NVAR,amb[1],Ja,1,&C[3][0],1);
1972  WAXPY(NVAR,amb[1],Jb,1,&C[3][0],1);
1973   
1974  Tau = T - a;
1975  WCOPY(NVAR,&C[3][0],1,Y,1);
1976  WSCAL(NVAR,pow(Tau,3),Y,1);
1977  for(j=2; j>=0; j--)
1978    WAXPY(NVAR,pow(Tau,(j-1)),&C[j][0],1,Y,1);
1979
1980} /* End of ros_Hermite3 */
1981
1982/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1983void ros_Hermite5( KPP_REAL a, KPP_REAL b, KPP_REAL T, KPP_REAL Ya[], 
1984                   KPP_REAL Yb[], KPP_REAL Ja[], KPP_REAL Jb[], KPP_REAL Ha[], 
1985                   KPP_REAL Hb[], KPP_REAL Y[] ) {
1986/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1987  Template for Hermite interpolation of order 5 on the interval [a,b]
1988 P = c(1) + c(2)*(x-a) + ... + c(6)*(x-a)^5
1989 P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb], P"[a,b] = [Ha,Hb]
1990~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1991
1992/*~~~> Local variables */
1993  KPP_REAL Tau, amb[5], C[6][NVAR];
1994  int i, j;
1995
1996  amb[0] = ((KPP_REAL)1.0)/(a-b);
1997  for(i=1; i<5; i++)
1998    amb[i] = amb[i-1]*amb[0];
1999
2000/* c(1) = ya; */
2001  WCOPY(NVAR,Ya,1,&C[0][0],1);
2002/* c(2) = ja; */
2003  WCOPY(NVAR,Ja,1,&C[1][0],1);
2004/* c(3) = ha/2; */
2005  WCOPY(NVAR,Ha,1,&C[2][0],1);
2006  WSCAL(NVAR,HALF,&C[2][0],1);
2007
2008/* c(4) = 10*amb(3)*ya - 10*amb(3)*yb - 6*amb(2)*ja - 4*amb(2)*jb 
2009          + 1.5*amb(1)*ha - 0.5*amb(1)*hb ; */
2010  WCOPY(NVAR,Ya,1,&C[3][0],1);
2011  WSCAL(NVAR,10.0*amb[2],&C[3][0],1);
2012  WAXPY(NVAR,-10.0*amb[2],Yb,1,&C[3][0],1);
2013  WAXPY(NVAR,-6.0*amb[1],Ja,1,&C[3][0],1);
2014  WAXPY(NVAR,-4.0*amb[1],Jb,1,&C[3][0],1);
2015  WAXPY(NVAR, 1.5*amb[0],Ha,1,&C[3][0],1);
2016  WAXPY(NVAR,-0.5*amb[0],Hb,1,&C[3][0],1);
2017
2018/* c(5) =   15*amb(4)*ya - 15*amb(4)*yb - 8.*amb(3)*ja - 7*amb(3)*jb
2019            + 1.5*amb(2)*ha - 1*amb(2)*hb ; */
2020  WCOPY(NVAR,Ya,1,&C[4][0],1);
2021  WSCAL(NVAR, 15.0*amb[3],&C[4][0],1);
2022  WAXPY(NVAR,-15.0*amb[3],Yb,1,&C[4][0],1);
2023  WAXPY(NVAR,-8.0*amb[2],Ja,1,&C[4][0],1);
2024  WAXPY(NVAR,-7.0*amb[2],Jb,1,&C[4][0],1);
2025  WAXPY(NVAR,1.5*amb[1],Ha,1,&C[4][0],1);
2026  WAXPY(NVAR,-amb[1],Hb,1,&C[4][0],1);
2027
2028/* c(6) =   6*amb(5)*ya - 6*amb(5)*yb - 3.*amb(4)*ja - 3.*amb(4)*jb
2029            + 0.5*amb(3)*ha -0.5*amb(3)*hb ; */
2030  WCOPY(NVAR,Ya,1,&C[5][0],1);
2031  WSCAL(NVAR, 6.0*amb[4],&C[5][0],1);
2032  WAXPY(NVAR,-6.0*amb[4],Yb,1,&C[5][0],1);
2033  WAXPY(NVAR,-3.0*amb[3],Ja,1,&C[5][0],1);
2034  WAXPY(NVAR,-3.0*amb[3],Jb,1,&C[5][0],1);
2035  WAXPY(NVAR, 0.5*amb[2],Ha,1,&C[5][0],1);
2036  WAXPY(NVAR,-0.5*amb[2],Hb,1,&C[5][0],1);
2037
2038  Tau = T - a;
2039  WCOPY(NVAR,&C[5][0],1,Y,1);
2040  for(j=4; j>=0; j--) {
2041    WSCAL(NVAR,Tau,Y,1);
2042    WAXPY(NVAR,ONE,&C[j][0],1,Y,1);
2043  }
2044
2045} /* End of ros_Hermite5 */
2046
2047/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2048void Ros2() {
2049/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2050 --- AN L-STABLE METHOD, 2 stages, order 2
2051~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2052
2053  KPP_REAL g;
2054
2055  g = (KPP_REAL)1.0 + ((KPP_REAL)1.0)/SQRT((KPP_REAL)2.0);
2056
2057  rosMethod = RS2;
2058/*~~~> Name of the method */
2059  strcpy(ros_Name, "ROS-2");
2060/*~~~> Number of stages */
2061  ros_S = 2;
2062
2063/*~~~> The coefficient matrices A and C are strictly lower triangular.
2064  The lower triangular (subdiagonal) elements are stored in row-wise order:
2065  A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc.
2066  The general mapping formula is:
2067      A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ]
2068      C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */
2069
2070  ros_A[0] = ((KPP_REAL)1.0)/g;
2071  ros_C[0] = ((KPP_REAL)-2.0)/g;
2072
2073/*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE)
2074       or does it re-use the function evaluation from stage i-1
2075       (ros_NewF[i]=FALSE) */
2076  ros_NewF[0] = TRUE;
2077  ros_NewF[1] = TRUE;
2078
2079/*~~~> M_i = Coefficients for new step solution */
2080  ros_M[0]= ((KPP_REAL)3.0)/((KPP_REAL)2.0*g);
2081  ros_M[1]= ((KPP_REAL)1.0)/((KPP_REAL)2.0*g);
2082
2083/* E_i = Coefficients for error estimator */
2084  ros_E[0] = ((KPP_REAL)1.0)/((KPP_REAL)2.0*g);
2085  ros_E[1] = ((KPP_REAL)1.0)/((KPP_REAL)2.0*g);
2086
2087/*~~~> ros_ELO = estimator of local order - the minimum between the
2088       main and the embedded scheme orders plus one */
2089  ros_ELO = (KPP_REAL)2.0;
2090
2091/*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
2092  ros_Alpha[0] = (KPP_REAL)0.0;
2093  ros_Alpha[1] = (KPP_REAL)1.0;
2094
2095/*~~~> Gamma_i = \sum_j  gamma_{i,j} */
2096  ros_Gamma[0] = g;
2097  ros_Gamma[1] =-g;
2098
2099} /* End of Ros2 */
2100
2101/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2102void Ros3() {
2103/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2104 --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
2105~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2106
2107  rosMethod = RS3;
2108/*~~~> Name of the method */
2109  strcpy(ros_Name, "ROS-3");
2110/*~~~> Number of stages */
2111  ros_S = 3;
2112
2113/*~~~> The coefficient matrices A and C are strictly lower triangular.
2114   The lower triangular (subdiagonal) elements are stored in row-wise order:
2115   A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc.
2116   The general mapping formula is:
2117       A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ]
2118       C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */
2119
2120  ros_A[0]= (KPP_REAL)1.0;
2121  ros_A[1]= (KPP_REAL)1.0;
2122  ros_A[2]= (KPP_REAL)0.0;
2123  ros_C[0] = (KPP_REAL)-0.10156171083877702091975600115545e01;
2124  ros_C[1] = (KPP_REAL) 0.40759956452537699824805835358067e01;
2125  ros_C[2] = (KPP_REAL) 0.92076794298330791242156818474003e01;
2126
2127/*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE)
2128       or does it re-use the function evaluation from stage i-1
2129       (ros_NewF[i]=FALSE) */
2130  ros_NewF[0] = TRUE;
2131  ros_NewF[1] = TRUE;
2132  ros_NewF[2] = FALSE;
2133/*~~~> M_i = Coefficients for new step solution */
2134  ros_M[0] = (KPP_REAL) 0.1e01;
2135  ros_M[1] = (KPP_REAL) 0.61697947043828245592553615689730e01;
2136  ros_M[2] = (KPP_REAL)-0.42772256543218573326238373806514;
2137/* E_i = Coefficients for error estimator */
2138  ros_E[0] = (KPP_REAL) 0.5;
2139  ros_E[1] = (KPP_REAL)-0.29079558716805469821718236208017e01;
2140  ros_E[2] = (KPP_REAL) 0.22354069897811569627360909276199;
2141
2142/*~~~> ros_ELO = estimator of local order - the minimum between the
2143       main and the embedded scheme orders plus 1 */
2144  ros_ELO = (KPP_REAL)3.0;
2145/*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
2146  ros_Alpha[0]= (KPP_REAL)0.0;
2147  ros_Alpha[1]= (KPP_REAL)0.43586652150845899941601945119356;
2148  ros_Alpha[2]= (KPP_REAL)0.43586652150845899941601945119356;
2149/*~~~> Gamma_i = \sum_j  gamma_{i,j} */
2150  ros_Gamma[0]= (KPP_REAL)0.43586652150845899941601945119356;
2151  ros_Gamma[1]= (KPP_REAL)0.24291996454816804366592249683314;
2152  ros_Gamma[2]= (KPP_REAL)0.21851380027664058511513169485832e01;
2153
2154} /* End of Ros3 */
2155
2156/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2157void Ros4() {
2158/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2159     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
2160     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
2161
2162      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2163      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2164      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2165      SPRINGER-VERLAG (1990)
2166~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2167
2168  rosMethod = RS4;
2169/*~~~> Name of the method */
2170  strcpy(ros_Name, "ROS-4");
2171/*~~~> Number of stages */
2172  ros_S = 4;
2173
2174/*~~~> The coefficient matrices A and C are strictly lower triangular.
2175   The lower triangular (subdiagonal) elements are stored in row-wise order:
2176   A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc.
2177   The general mapping formula is:
2178       A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ]
2179       C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */
2180
2181  ros_A[0] = (KPP_REAL)0.2000000000000000e01;
2182  ros_A[1] = (KPP_REAL)0.1867943637803922e01;
2183  ros_A[2] = (KPP_REAL)0.2344449711399156;
2184  ros_A[3] = ros_A[1];
2185  ros_A[4] = ros_A[2];
2186  ros_A[5] = (KPP_REAL)0.0;
2187
2188  ros_C[0] = (KPP_REAL)-0.7137615036412310e01;
2189  ros_C[1] = (KPP_REAL) 0.2580708087951457e01;
2190  ros_C[2] = (KPP_REAL) 0.6515950076447975;
2191  ros_C[3] = (KPP_REAL)-0.2137148994382534e01;
2192  ros_C[4] = (KPP_REAL)-0.3214669691237626;
2193  ros_C[5] = (KPP_REAL)-0.6949742501781779;
2194
2195/*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE)
2196       or does it re-use the function evaluation from stage i-1
2197       (ros_NewF[i]=FALSE) */
2198  ros_NewF[0] = TRUE;
2199  ros_NewF[1] = TRUE;
2200  ros_NewF[2] = TRUE;
2201  ros_NewF[3] = FALSE;
2202/*~~~> M_i = Coefficients for new step solution */
2203  ros_M[0] = (KPP_REAL)0.2255570073418735e01;
2204  ros_M[1] = (KPP_REAL)0.2870493262186792;
2205  ros_M[2] = (KPP_REAL)0.4353179431840180;
2206  ros_M[3] = (KPP_REAL)0.1093502252409163e01;
2207/*~~~> E_i  = Coefficients for error estimator */
2208  ros_E[0] = (KPP_REAL)-0.2815431932141155;
2209  ros_E[1] = (KPP_REAL)-0.7276199124938920e-01;
2210  ros_E[2] = (KPP_REAL)-0.1082196201495311;
2211  ros_E[3] = (KPP_REAL)-0.1093502252409163e01;
2212
2213/*~~~> ros_ELO  = estimator of local order - the minimum between the
2214       main and the embedded scheme orders plus 1 */
2215  ros_ELO = (KPP_REAL)4.0;
2216/*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
2217  ros_Alpha[0] = (KPP_REAL)0.0;
2218  ros_Alpha[1] = (KPP_REAL)0.1145640000000000e01;
2219  ros_Alpha[2] = (KPP_REAL)0.6552168638155900;
2220  ros_Alpha[3] = ros_Alpha[2];
2221/*~~~> Gamma_i = \sum_j  gamma_{i,j} */ 
2222  ros_Gamma[0] = (KPP_REAL) 0.5728200000000000;
2223  ros_Gamma[1] = (KPP_REAL)-0.1769193891319233e01;
2224  ros_Gamma[2] = (KPP_REAL) 0.7592633437920482;
2225  ros_Gamma[3] = (KPP_REAL)-0.1049021087100450;
2226
2227} /* End of Ros4 */
2228
2229/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2230void Rodas3() {
2231/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2232 --- A STIFFLY-STABLE METHOD, 4 stages, order 3
2233~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2234
2235  rosMethod = RD3;
2236/*~~~> Name of the method */
2237  strcpy(ros_Name, "RODAS-3");
2238/*~~~> Number of stages */
2239  ros_S = 4;
2240
2241/*~~~> The coefficient matrices A and C are strictly lower triangular.
2242   The lower triangular (subdiagonal) elements are stored in row-wise order:
2243   A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc.
2244   The general mapping formula is:
2245       A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ]
2246       C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */
2247
2248  ros_A[0] = (KPP_REAL)0.0;
2249  ros_A[1] = (KPP_REAL)2.0;
2250  ros_A[2] = (KPP_REAL)0.0;
2251  ros_A[3] = (KPP_REAL)2.0;
2252  ros_A[4] = (KPP_REAL)0.0;
2253  ros_A[5] = (KPP_REAL)1.0;
2254
2255  ros_C[0] = (KPP_REAL) 4.0;
2256  ros_C[1] = (KPP_REAL) 1.0;
2257  ros_C[2] = (KPP_REAL)-1.0;
2258  ros_C[3] = (KPP_REAL) 1.0;
2259  ros_C[4] = (KPP_REAL)-1.0;
2260  ros_C[5] = -(((KPP_REAL)8.0)/((KPP_REAL)3.0));
2261
2262/*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE)
2263       or does it re-use the function evaluation from stage i-1
2264       (ros_NewF[i]=FALSE) */
2265  ros_NewF[0] = TRUE;
2266  ros_NewF[1] = FALSE;
2267  ros_NewF[2] = TRUE;
2268  ros_NewF[3] = TRUE;
2269/*~~~> M_i = Coefficients for new step solution */
2270  ros_M[0] = (KPP_REAL)2.0;
2271  ros_M[1] = (KPP_REAL)0.0;
2272  ros_M[2] = (KPP_REAL)1.0;
2273  ros_M[3] = (KPP_REAL)1.0;
2274/*~~~> E_i  = Coefficients for error estimator */
2275  ros_E[0] = (KPP_REAL)0.0;
2276  ros_E[1] = (KPP_REAL)0.0;
2277  ros_E[2] = (KPP_REAL)0.0;
2278  ros_E[3] = (KPP_REAL)1.0;
2279
2280/*~~~> ros_ELO  = estimator of local order - the minimum between the
2281    main and the embedded scheme orders plus 1 */
2282  ros_ELO  = (KPP_REAL)3.0;
2283/*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
2284  ros_Alpha[0] = (KPP_REAL)0.0;
2285  ros_Alpha[1] = (KPP_REAL)0.0;
2286  ros_Alpha[2] = (KPP_REAL)1.0;
2287  ros_Alpha[3] = (KPP_REAL)1.0;
2288/*~~~> Gamma_i = \sum_j  gamma_{i,j} */
2289  ros_Gamma[0] = (KPP_REAL)0.5;
2290  ros_Gamma[1] = (KPP_REAL)1.5;
2291  ros_Gamma[2] = (KPP_REAL)0.0;
2292  ros_Gamma[3] = (KPP_REAL)0.0;
2293
2294} /* End of Rodas3 */
2295
2296/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2297void Rodas4() {
2298/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2299     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
2300
2301      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2302      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2303      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2304      SPRINGER-VERLAG (1996)
2305~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2306
2307  rosMethod = RD4;
2308/*~~~> Name of the method */
2309  strcpy(ros_Name, "RODAS-4");
2310/*~~~> Number of stages */
2311  ros_S = 6;
2312
2313/*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */
2314  ros_Alpha[0] = (KPP_REAL)0.000;
2315  ros_Alpha[1] = (KPP_REAL)0.386;
2316  ros_Alpha[2] = (KPP_REAL)0.210;
2317  ros_Alpha[3] = (KPP_REAL)0.630;
2318  ros_Alpha[4] = (KPP_REAL)1.000;
2319  ros_Alpha[5] = (KPP_REAL)1.000;
2320
2321/*~~~> Gamma_i = \sum_j  gamma_{i,j} */
2322  ros_Gamma[0] = (KPP_REAL) 0.2500000000000000;
2323  ros_Gamma[1] = (KPP_REAL)-0.1043000000000000;
2324  ros_Gamma[2] = (KPP_REAL) 0.1035000000000000;
2325  ros_Gamma[3] = (KPP_REAL)-0.3620000000000023e-01;
2326  ros_Gamma[4] = (KPP_REAL) 0.0;
2327  ros_Gamma[5] = (KPP_REAL) 0.0;
2328
2329/*~~~> The coefficient matrices A and C are strictly lower triangular.
2330   The lower triangular (subdiagonal) elements are stored in row-wise order:
2331   A[0][1] = ros_A[0], A[0][2]=ros_A[1], A[1][2]=ros_A[2], etc.
2332   The general mapping formula is:  A[i][j] = ros_A[ (i-1)*(i-2)/2 + j ]
2333                  C[i][j] = ros_C[ (i-1)*(i-2)/2 + j ] */
2334
2335  ros_A[0] = (KPP_REAL) 0.1544000000000000e01;
2336  ros_A[1] = (KPP_REAL) 0.9466785280815826;
2337  ros_A[2] = (KPP_REAL) 0.2557011698983284;
2338  ros_A[3] = (KPP_REAL) 0.3314825187068521e01;
2339  ros_A[4] = (KPP_REAL) 0.2896124015972201e01;
2340  ros_A[5] = (KPP_REAL) 0.9986419139977817;
2341  ros_A[6] = (KPP_REAL) 0.1221224509226641e01;
2342  ros_A[7] = (KPP_REAL) 0.6019134481288629e01;
2343  ros_A[8] = (KPP_REAL) 0.1253708332932087e02;
2344  ros_A[9] = (KPP_REAL)-0.6878860361058950;
2345  ros_A[10] = ros_A[6];
2346  ros_A[11] = ros_A[7];
2347  ros_A[12] = ros_A[8];
2348  ros_A[13] = ros_A[9];
2349  ros_A[14] = (KPP_REAL)1.0;
2350 
2351  ros_C[0]  = (KPP_REAL)-0.5668800000000000e01;
2352  ros_C[1]  = (KPP_REAL)-0.2430093356833875e01;
2353  ros_C[2]  = (KPP_REAL)-0.2063599157091915;
2354  ros_C[3]  = (KPP_REAL)-0.1073529058151375;
2355  ros_C[4]  = (KPP_REAL)-0.9594562251023355e01;
2356  ros_C[5]  = (KPP_REAL)-0.2047028614809616e02;
2357  ros_C[6]  = (KPP_REAL) 0.7496443313967647e01;
2358  ros_C[7]  = (KPP_REAL)-0.1024680431464352e02;
2359  ros_C[8]  = (KPP_REAL)-0.3399990352819905e02;
2360  ros_C[9]  = (KPP_REAL) 0.1170890893206160e02;
2361  ros_C[10] = (KPP_REAL) 0.8083246795921522e01;
2362  ros_C[11] = (KPP_REAL)-0.7981132988064893e01;
2363  ros_C[12] = (KPP_REAL)-0.3152159432874371e02;
2364  ros_C[13] = (KPP_REAL) 0.1631930543123136e02;
2365  ros_C[14] = (KPP_REAL)-0.6058818238834054e01;
2366
2367/*~~~> M_i = Coefficients for new step solution */
2368  ros_M[0] = ros_A[6];
2369  ros_M[1] = ros_A[7];
2370  ros_M[2] = ros_A[8];
2371  ros_M[3] = ros_A[9];
2372  ros_M[4] = (KPP_REAL)1.0;
2373  ros_M[5] = (KPP_REAL)1.0;
2374
2375/*~~~> E_i  = Coefficients for error estimator */
2376  ros_E[0] = (KPP_REAL)0.0;
2377  ros_E[1] = (KPP_REAL)0.0;
2378  ros_E[2] = (KPP_REAL)0.0;
2379  ros_E[3] = (KPP_REAL)0.0;
2380  ros_E[4] = (KPP_REAL)0.0;
2381  ros_E[5] = (KPP_REAL)1.0;
2382
2383/*~~~> Does the stage i require a new function evaluation (ros_NewF[i]=TRUE)
2384       or does it re-use the function evaluation from stage i-1
2385       (ros_NewF[i]=FALSE) */
2386  ros_NewF[0] = TRUE;
2387  ros_NewF[1] = TRUE;
2388  ros_NewF[2] = TRUE;
2389  ros_NewF[3] = TRUE;
2390  ros_NewF[4] = TRUE;
2391  ros_NewF[5] = TRUE;
2392
2393/*~~~> ros_ELO  = estimator of local order - the minimum between the
2394        main and the embedded scheme orders plus 1 */
2395  ros_ELO = (KPP_REAL)4.0;
2396
2397} /* End of Rodas4 */
2398
2399/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2400void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ) {
2401/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2402  Template for the ODE function call.
2403  Updates the rate coefficients (and possibly the fixed species) at each call
2404~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2405
2406/*~~~> Local variables */
2407  KPP_REAL Told;
2408
2409  Told = TIME;
2410  TIME = T;
2411  Update_SUN();
2412  Update_RCONST();
2413  Fun( Y, FIX, RCONST, Ydot );
2414  TIME = Told;
2415
2416} /* End of FunTemplate */
2417
2418/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2419void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] ) {
2420/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2421  Template for the ODE Jacobian call.
2422  Updates the rate coefficients (and possibly the fixed species) at each call
2423~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2424
2425/*~~~> Local variables */
2426  KPP_REAL Told;
2427#ifdef FULL_ALGEBRA
2428  KPP_REAL JV[LU_NONZERO];
2429  int i, j;
2430#endif
2431
2432  Told = TIME;
2433  TIME = T;
2434  Update_SUN();
2435  Update_RCONST();
2436#ifdef FULL_ALGEBRA
2437  Jac_SP(Y, FIX, RCONST, JV);
2438  for(j=0; j<NVAR; j++) {
2439    for(i=0; i<NVAR; i++)
2440      Jcb[i][j] = (KPP_REAL)0.0;
2441  }
2442  for(i=0; i<LU_NONZERO; i++)
2443    Jcb[LU_ICOL[i]][LU_IROW[i]] = JV[i];
2444#else
2445  Jac_SP( Y, FIX, RCONST, Jcb );
2446#endif
2447  TIME = Told;
2448} /* End of JacTemplate */
2449
2450/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2451void HessTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Hes[] ) {
2452/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2453  Template for the ODE Hessian call.
2454  Updates the rate coefficients (and possibly the fixed species) at each call
2455~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2456
2457/*~~~> Local variables */
2458  KPP_REAL Told;
2459
2460  Told = TIME;
2461  TIME = T;
2462  Update_SUN();
2463  Update_RCONST();
2464  Hessian( Y, FIX, RCONST, Hes );
2465  TIME = Told;
2466
2467} /* End of HessTemplate */
2468
2469/* End of INTEGRATE function                                                 */
2470/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |