Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock.f @ 3698

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

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/rosenbrock.f @ 3698

Download in other formats: