Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

runge_kutta.c @ 4218

Last change on this file since 4218 was 2696, checked in by kanani, 7 years ago
Merge of branch palm4u into trunk
File size: 75.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: