Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module mod_dopri5
21 : use const_def, only: dp
22 : use math_lib
23 :
24 : implicit none
25 :
26 : contains
27 :
28 11 : subroutine do_dopri5(
29 : & n,fcn,x,y,xend,h,max_step_size,max_steps,
30 11 : & rtol,atol,itol,solout,iout,work,lwork,iwork,liwork,
31 : & lrpar,rpar,lipar,ipar,lout,idid)
32 : ! *** *** *** *** *** *** *** *** *** *** *** *** ***
33 : ! declarations
34 : ! *** *** *** *** *** *** *** *** *** *** *** *** ***
35 : integer, intent(in) :: n ! the dimension of the system
36 : interface
37 : #include "num_fcn.dek"
38 : #include "num_solout.dek"
39 : end interface
40 : real(dp), intent(inout) :: x
41 : real(dp), intent(inout), pointer :: y(:) ! (n)
42 : real(dp), intent(in) :: xend, h, max_step_size
43 : real(dp), intent(in) :: rtol(*)
44 : real(dp), intent(in) :: atol(*)
45 : integer, intent(in) :: itol, iout, liwork, lwork, max_steps
46 : integer, intent(inout) :: iwork(liwork)
47 : real(dp), intent(inout) :: work(lwork)
48 : integer, intent(in) :: lrpar, lipar
49 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
50 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
51 : integer, intent(in) :: lout
52 : integer, intent(out) :: idid
53 : integer :: nfcn, nstep, naccpt, nrejct, i, icomp, ieco
54 : integer :: iek1, iek2, iek3, iek4, iek5, iek6
55 : integer :: iey1, ieys, istore, nmax, meth, nstiff, nrdens
56 11 : real(dp) :: beta, fac1, fac2, hmax, safe, uround
57 :
58 : logical :: arret
59 : ! *** *** *** *** *** *** ***
60 : ! setting the parameters
61 : ! *** *** *** *** *** *** ***
62 11 : nfcn=0
63 11 : nstep=0
64 11 : naccpt=0
65 11 : nrejct=0
66 11 : arret=.false.
67 : ! -------- nmax , the maximal number of steps -----
68 11 : if (max_steps == 0 ) then
69 1 : nmax=100000
70 : else
71 10 : nmax=max_steps
72 10 : if (nmax <= 0 ) then
73 0 : if (lout > 0) write(lout,*) ' wrong input max_steps=',max_steps
74 : arret=.true.
75 : end if
76 : end if
77 : ! -------- meth coefficients of the method
78 11 : if (iwork(2) == 0 ) then
79 11 : meth=1
80 : else
81 0 : meth=iwork(2)
82 0 : if (meth <= 0.or.meth >= 4 ) then
83 0 : if (lout > 0) write(lout,*) ' curious input iwork(2)=',iwork(2)
84 : arret=.true.
85 : end if
86 : end if
87 : ! -------- nstiff parameter for stiffness detection
88 11 : nstiff=iwork(4)
89 11 : if (nstiff == 0) nstiff=1000
90 11 : if (nstiff < 0) nstiff=nmax+10
91 : ! -------- nrdens number of dense output components
92 11 : nrdens=iwork(5)
93 11 : if (nrdens < 0.or.nrdens > n ) then
94 0 : if (lout > 0) write(lout,*) ' curious input iwork(5)=',iwork(5)
95 : arret=.true.
96 : else
97 11 : if (nrdens > 0.and.iout < 2 ) then
98 0 : if (lout > 0) write(lout,*) ' warning: put iout=2 for dense output '
99 : end if
100 11 : if (nrdens == n) then
101 3 : do i=1,nrdens
102 3 : iwork(20+i)=i
103 : end do
104 : end if
105 : end if
106 : ! -------- uround smallest number satisfying 1.d0+uround>1.d0
107 11 : if (work(1) == 0.d0 ) then
108 11 : uround=2.3d-16
109 : else
110 0 : uround=work(1)
111 0 : if (uround <= 1.d-35.or.uround >= 1.d0 ) then
112 0 : if (lout > 0) write(lout,*) ' which machine do you have? your uround was:',work(1)
113 : arret=.true.
114 : end if
115 : end if
116 : ! ------- safety factor -------------
117 11 : if (work(2) == 0.d0 ) then
118 11 : safe=0.9d0
119 : else
120 0 : safe=work(2)
121 0 : if (safe >= 1.d0.or.safe <= 1.d-4 ) then
122 0 : if (lout > 0) write(lout,*) ' curious input for safety factor work(2)=',work(2)
123 : arret=.true.
124 : end if
125 : end if
126 : ! ------- fac1,fac2 parameters for step size selection
127 11 : if (work(3) == 0.d0 ) then
128 11 : fac1=0.2d0
129 : else
130 0 : fac1=work(3)
131 : end if
132 11 : if (work(4) == 0.d0 ) then
133 11 : fac2=10.d0
134 : else
135 0 : fac2=work(4)
136 : end if
137 : ! --------- beta for step control stabilization -----------
138 11 : if (work(5) == 0.d0 ) then
139 11 : beta=0.04d0
140 : else
141 0 : if (work(5) < 0.d0 ) then
142 0 : beta=0.d0
143 : else
144 0 : beta=work(5)
145 0 : if (beta > 0.2d0 ) then
146 0 : if (lout > 0) write(lout,*) ' curious input for beta: work(5)=',work(5)
147 : arret=.true.
148 : end if
149 : end if
150 : end if
151 : ! -------- maximal step size
152 11 : if (max_step_size == 0.d0 ) then
153 10 : hmax=xend-x
154 : else
155 1 : hmax=max_step_size
156 : end if
157 : ! ------- prepare the entry-points for the arrays in work -----
158 11 : iey1=21
159 11 : iek1=iey1+n
160 11 : iek2=iek1+n
161 11 : iek3=iek2+n
162 11 : iek4=iek3+n
163 11 : iek5=iek4+n
164 11 : iek6=iek5+n
165 11 : ieys=iek6+n
166 11 : ieco=ieys+n
167 : ! ------ total storage requirement -----------
168 11 : istore=ieys+(3+5*nrdens)-1
169 11 : if (istore > lwork ) then
170 0 : if (lout > 0) write(lout,*) ' insufficient storage for work, min. lwork=',istore
171 : arret=.true.
172 : end if
173 11 : icomp=21
174 11 : istore=icomp+nrdens-1
175 11 : if (istore > liwork ) then
176 0 : if (lout > 0) write(lout,*) ' insufficient storage for iwork, min. liwork=',istore
177 : arret=.true.
178 : end if
179 : ! ------ when a fail has occurred, we return with idid=-1
180 11 : if (arret) then
181 0 : idid=-1
182 0 : return
183 : end if
184 : ! -------- call to core integrator ------------
185 : call dopcor(n,fcn,x,y,xend,hmax,h,rtol,atol,itol,lout,
186 : & solout,iout,idid,nmax,uround,meth,nstiff,safe,beta,fac1,fac2,
187 : & work(iey1),work(iek1),work(iek2),work(iek3),work(iek4),
188 : & work(iek5),work(iek6),work(ieys),work(ieco),iwork(icomp),
189 11 : & nrdens,lrpar,rpar,lipar,ipar,nfcn,nstep,naccpt,nrejct)
190 11 : work(7)=h
191 11 : iwork(17)=nfcn
192 11 : iwork(18)=nstep
193 11 : iwork(19)=naccpt
194 11 : iwork(20)=nrejct
195 : ! ----------- return -----------
196 11 : return
197 : end subroutine do_dopri5
198 :
199 :
200 : !
201 : ! ----- ... and here is the core integrator ----------
202 : !
203 11 : subroutine dopcor(n,fcn,x,y,xend,hmax,h,rtol,atol,itol,lout,
204 : & solout,iout,idid,nmax,uround,meth,nstiff,safe,beta,fac1,fac2,
205 11 : & y1,k1,k2,k3,k4,k5,k6,ysti,rwork,icomp,nrd,lrpar,rpar,lipar,ipar,
206 : & nfcn,nstep,naccpt,nrejct)
207 : ! ----------------------------------------------------------
208 : ! core integrator for dopri5
209 : ! parameters same as in dopri5 with workspace added
210 : ! ----------------------------------------------------------
211 : ! declarations
212 : ! ----------------------------------------------------------
213 : integer :: n, nrd, itol, lout, iout, idid, nmax, meth
214 11 : real(dp) :: x, xold, xend, hmax, h, uround, safe
215 : real(dp) :: beta, fac1, fac2
216 : real(dp) :: k1(n),k2(n),k3(n),k4(n),k5(n),k6(n)
217 : real(dp) :: y(n),y1(n),ysti(n),atol(*),rtol(*)
218 22 : integer :: icomp(nrd), iwork(nrd+1)
219 : integer :: lrpar, lipar, nfcn, nstep, naccpt, nrejct, nstiff
220 : logical :: reject,last
221 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
222 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
223 11 : real(dp) :: c2,c3,c4,c5
224 11 : real(dp) :: e1,e3,e4,e5,e6,e7
225 11 : real(dp) :: a21,a31,a32,a41,a42,a43,a51,a52,a53,a54
226 11 : real(dp) :: a61,a62,a63,a64,a65,a71,a73,a74,a75,a76
227 11 : real(dp) :: d1,d3,d4,d5,d6,d7
228 11 : real(dp) :: atoli, expo1, facc1, facc2, facold, hlamb, hout, nonsti, posneg, bspl, rtoli
229 : integer :: i, iasti, ierr, iord, irtrn, j
230 : !common /condo5/xold,hout
231 11 : real(dp) :: err, fac, fac11, hnew, sk, stden, stnum, xph, yd0, ydiff
232 :
233 : interface
234 : #include "num_fcn.dek"
235 : #include "num_solout.dek"
236 : end interface
237 :
238 : real(dp), target :: rwork(3+5*nrd)
239 11 : real(dp), pointer :: cont(:)
240 11 : cont => rwork(3:3+5*nrd)
241 :
242 :
243 : ! *** *** *** *** *** *** ***
244 : ! initializations
245 : ! *** *** *** *** *** *** ***
246 11 : if (meth == 1) call cdopri(c2,c3,c4,c5,e1,e3,e4,e5,e6,e7,
247 : & a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
248 : & a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,
249 11 : & d1,d3,d4,d5,d6,d7)
250 11 : facold=1.d-4
251 11 : expo1=0.2d0-beta*0.75d0
252 11 : facc1=1.d0/fac1
253 11 : facc2=1.d0/fac2
254 11 : posneg=sign(1.d0,xend-x)
255 11 : nonsti=0
256 : ! --- initial preparations
257 11 : atoli=atol(1)
258 11 : rtoli=rtol(1)
259 11 : last=.false.
260 11 : hlamb=0.d0
261 11 : iasti=0
262 11 : call fcn(n,x,h,y,k1,lrpar,rpar,lipar,ipar,ierr)
263 11 : if (ierr /= 0) then; idid=-5; return; end if
264 11 : hmax=abs(hmax)
265 11 : iord=5
266 11 : if (h == 0.d0) h=hinit(n,fcn,x,y,xend,posneg,k1,k2,k3,iord,
267 0 : & hmax,atol,rtol,itol,lrpar,rpar,lipar,ipar,ierr)
268 11 : if (ierr /= 0) then; idid=-5; return; end if
269 11 : nfcn=nfcn+2
270 11 : reject=.false.
271 11 : xold=x
272 11 : if (iout /= 0) then
273 1 : irtrn=1
274 1 : hout=h
275 1 : rwork(1) = xold
276 1 : rwork(2) = hout
277 1 : iwork(1) = nrd
278 3 : iwork(2:nrd+1) = icomp(1:nrd)
279 :
280 1 : if (iout >= 2) then
281 3 : do j=1,nrd
282 2 : i=icomp(j)
283 2 : cont(j)=y(i)
284 2 : cont(nrd+j)=0
285 2 : cont(2*nrd+j)=0
286 3 : cont(3*nrd+j)=0
287 : end do
288 : end if
289 :
290 1 : call solout(naccpt+1,xold,x,n,y,rwork,iwork,contd5,lrpar,rpar,lipar,ipar,irtrn)
291 1 : if (irtrn < 0) GOTO 79
292 : else
293 10 : irtrn=0
294 : end if
295 : ! --- basic integration step
296 : 1 continue
297 2266 : if (nstep > nmax) GOTO 78
298 2266 : if (0.1d0*abs(h) <= abs(x)*uround) GOTO 77
299 2266 : if ((x+1.01d0*h-xend)*posneg > 0.d0) then
300 11 : h=xend-x
301 11 : last=.true.
302 : end if
303 2266 : nstep=nstep+1
304 : ! --- the first 6 stages
305 2266 : if (irtrn >= 2) then
306 0 : call fcn(n,x,h,y,k1,lrpar,rpar,lipar,ipar,ierr)
307 0 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
308 : end if
309 5770 : do i=1,n
310 5770 : y1(i)=y(i)+h*a21*k1(i)
311 : end do
312 2266 : call fcn(n,x+c2*h,h,y1,k2,lrpar,rpar,lipar,ipar,ierr)
313 2266 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
314 5770 : do i=1,n
315 5770 : y1(i)=y(i)+h*(a31*k1(i)+a32*k2(i))
316 : end do
317 2266 : call fcn(n,x+c3*h,h,y1,k3,lrpar,rpar,lipar,ipar,ierr)
318 2266 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
319 5770 : do i=1,n
320 5770 : y1(i)=y(i)+h*(a41*k1(i)+a42*k2(i)+a43*k3(i))
321 : end do
322 2266 : call fcn(n,x+c4*h,h,y1,k4,lrpar,rpar,lipar,ipar,ierr)
323 2266 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
324 5770 : do i=1,n
325 5770 : y1(i)=y(i)+h*(a51*k1(i)+a52*k2(i)+a53*k3(i)+a54*k4(i))
326 : end do
327 2266 : call fcn(n,x+c5*h,h,y1,k5,lrpar,rpar,lipar,ipar,ierr)
328 2266 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
329 5770 : do i=1,n
330 5770 : ysti(i)=y(i)+h*(a61*k1(i)+a62*k2(i)+a63*k3(i)+a64*k4(i)+a65*k5(i))
331 : end do
332 2266 : xph=x+h
333 2266 : call fcn(n,xph,h,ysti,k6,lrpar,rpar,lipar,ipar,ierr)
334 2266 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
335 5770 : do i=1,n
336 5770 : y1(i)=y(i)+h*(a71*k1(i)+a73*k3(i)+a74*k4(i)+a75*k5(i)+a76*k6(i))
337 : end do
338 2266 : call fcn(n,xph,h,y1,k2,lrpar,rpar,lipar,ipar,ierr)
339 2266 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
340 2266 : if (iout >= 2) then
341 3714 : do j=1,nrd
342 2476 : i=icomp(j)
343 : cont(4*nrd+j)=h*(d1*k1(i)+d3*k3(i)+d4*k4(i)+d5*k5(i)
344 3714 : & +d6*k6(i)+d7*k2(i))
345 : end do
346 : end if
347 5770 : do i=1,n
348 5770 : k4(i)=(e1*k1(i)+e3*k3(i)+e4*k4(i)+e5*k5(i)+e6*k6(i)+e7*k2(i))*h
349 : end do
350 2266 : nfcn=nfcn+6
351 : ! --- error estimation
352 2266 : err=0.d0
353 2266 : if (itol == 0) then
354 3714 : do i=1,n
355 2476 : sk=atoli+rtoli*max(abs(y(i)),abs(y1(i)))
356 3714 : err=err+pow2(k4(i)/sk)
357 : end do
358 : else
359 2056 : do i=1,n
360 1028 : sk=atol(i)+rtol(i)*max(abs(y(i)),abs(y1(i)))
361 2056 : err=err+pow2(k4(i)/sk)
362 : end do
363 : end if
364 2266 : err=sqrt(err/n)
365 : ! --- computation of hnew
366 2266 : fac11=pow(err,expo1)
367 : ! --- lund-stabilization
368 2266 : fac=fac11/pow(facold,beta)
369 : ! --- we require fac1 <= hnew/h <= fac2
370 2266 : fac=max(facc2,min(facc1,fac/safe))
371 2266 : hnew=h/fac
372 2266 : if (err <= 1.d0 ) then
373 : ! --- step is accepted
374 2201 : facold=max(err,1.0d-4)
375 2201 : naccpt=naccpt+1
376 : ! ------- stiffness detection
377 2201 : if (mod(naccpt,nstiff) == 0.or.iasti > 0) then
378 : stnum=0.d0
379 : stden=0.d0
380 3639 : do i=1,n
381 2426 : stnum=stnum+pow2(k2(i)-k6(i))
382 3639 : stden=stden+pow2(y1(i)-ysti(i))
383 : end do
384 1213 : if (stden > 0.d0) hlamb=h*sqrt(stnum/stden)
385 1213 : if (hlamb > 3.25d0) then
386 1076 : nonsti=0
387 1076 : iasti=iasti+1
388 1076 : if (iasti == 15) then
389 3 : if (lout > 0) write (lout,*) ' the problem seems to become stiff at x = ',x
390 3 : if (lout < 0) GOTO 76
391 : end if
392 : else
393 137 : nonsti=nonsti+1
394 137 : if (nonsti == 6) iasti=0
395 : end if
396 : end if
397 2201 : if (iout >= 2) then
398 3639 : do j=1,nrd
399 2426 : i=icomp(j)
400 2426 : yd0=y(i)
401 2426 : ydiff=y1(i)-yd0
402 2426 : bspl=h*k1(i)-ydiff
403 2426 : cont(j)=y(i)
404 2426 : cont(nrd+j)=ydiff
405 2426 : cont(2*nrd+j)=bspl
406 3639 : cont(3*nrd+j)=-h*k2(i)+ydiff-bspl
407 : end do
408 : end if
409 5615 : do i=1,n
410 3414 : k1(i)=k2(i)
411 5615 : y(i)=y1(i)
412 : end do
413 2201 : xold=x
414 2201 : x=xph
415 2201 : if (iout /= 0) then
416 : irtrn=1
417 1213 : hout=h
418 1213 : rwork(1) = xold
419 1213 : rwork(2) = hout
420 1213 : iwork(1) = nrd
421 3639 : iwork(2:nrd+1) = icomp(1:nrd)
422 1213 : call solout(naccpt+1,xold,x,n,y,rwork,iwork,contd5,lrpar,rpar,lipar,ipar,irtrn)
423 1213 : if (irtrn < 0) GOTO 79
424 : end if
425 : ! ------- normal exit
426 2201 : if (last) then
427 11 : h=hnew
428 11 : idid=1
429 11 : return
430 : end if
431 2190 : if (abs(hnew) > hmax)hnew=posneg*hmax
432 2190 : if (reject)hnew=posneg*min(abs(hnew),abs(h))
433 : reject=.false.
434 : else
435 : ! --- step is rejected
436 65 : hnew=h/min(facc1,fac11/safe)
437 65 : reject=.true.
438 65 : if (naccpt >= 1)nrejct=nrejct+1
439 : last=.false.
440 : end if
441 2255 : h=hnew
442 2255 : GOTO 1
443 : ! --- fail exit
444 : 76 continue
445 0 : idid=-4
446 0 : return
447 : 77 continue
448 0 : if (lout > 0) write(lout,979)x
449 0 : if (lout > 0) write(lout,*)' step size too small, h=',h
450 0 : idid=-3
451 0 : return
452 : 78 continue
453 0 : if (lout > 0) write(lout,979)x
454 0 : if (lout > 0) write(lout,*) ' more than nmax =',nmax,'steps are needed'
455 0 : idid=-2
456 0 : return
457 : 79 continue
458 : !if (lout > 0) write(lout,979)x
459 : 979 format(' exit of dopri5 at x=',e18.4)
460 0 : idid=2
461 0 : return
462 11 : end subroutine dopcor
463 : !
464 0 : real(dp) function hinit(n,fcn,x,y,xend,posneg,f0,f1,y1,iord,
465 : & hmax,atol,rtol,itol,lrpar,rpar,lipar,ipar,ierr)
466 : ! ----------------------------------------------------------
467 : ! ---- computation of an initial step size guess
468 : ! ----------------------------------------------------------
469 : integer, intent(in) :: n
470 : real(dp) :: x
471 : dimension y(n),y1(n),f0(n),f1(n),atol(*),rtol(*)
472 : integer, intent(in) :: lrpar, lipar
473 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
474 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
475 0 : real(dp) :: y, xend, posneg, f0, f1, y1, hmax, atol, rtol, atoli
476 0 : real(dp) :: der2, der12, dnf, dny, h, h1, rtoli, sk
477 : integer :: i, iord, itol, ierr, idid
478 :
479 : interface
480 : #include "num_fcn.dek"
481 : end interface
482 :
483 : ! ---- compute a first guess for explicit euler as
484 : ! ---- h = 0.01 * norm (y0) / norm (f0)
485 : ! ---- the increment for explicit euler is small
486 : ! ---- compared to the solution
487 0 : dnf=0.0d0
488 0 : dny=0.0d0
489 0 : atoli=atol(1)
490 0 : rtoli=rtol(1)
491 0 : if (itol == 0) then
492 0 : do i=1,n
493 0 : sk=atoli+rtoli*abs(y(i))
494 0 : dnf=dnf+pow2(f0(i)/sk)
495 0 : dny=dny+pow2(y(i)/sk)
496 : end do
497 : else
498 0 : do i=1,n
499 0 : sk=atol(i)+rtol(i)*abs(y(i))
500 0 : dnf=dnf+pow2(f0(i)/sk)
501 0 : dny=dny+pow2(y(i)/sk)
502 : end do
503 : end if
504 0 : if (dnf <= 1.d-10.or.dny <= 1.d-10) then
505 0 : h=1.0d-6
506 : else
507 0 : h=sqrt(dny/dnf)*0.01d0
508 : end if
509 0 : h=min(h,hmax)
510 0 : h=sign(h,posneg)
511 : ! ---- perform an explicit euler step
512 0 : do i=1,n
513 0 : y1(i)=y(i)+h*f0(i)
514 : end do
515 0 : call fcn(n,x+h,h,y1,f1,lrpar,rpar,lipar,ipar,ierr)
516 0 : if (ierr /= 0) then; idid=-5; return; end if
517 : ! ---- estimate the second derivative of the solution
518 0 : der2=0.0d0
519 0 : if (itol == 0) then
520 0 : do i=1,n
521 0 : sk=atoli+rtoli*abs(y(i))
522 0 : der2=der2+pow2((f1(i)-f0(i))/sk)
523 : end do
524 : else
525 0 : do i=1,n
526 0 : sk=atol(i)+rtol(i)*abs(y(i))
527 0 : der2=der2+pow2((f1(i)-f0(i))/sk)
528 : end do
529 : end if
530 0 : der2=sqrt(der2)/h
531 : ! ---- step size is computed such that
532 : ! ---- h**iord * max ( norm (f0), norm (der2)) = 0.01
533 0 : der12=max(abs(der2),sqrt(dnf))
534 0 : if (der12 <= 1.d-15) then
535 0 : h1=max(1.0d-6,abs(h)*1.0d-3)
536 : else
537 0 : h1=pow(0.01d0/der12,1.d0/iord)
538 : end if
539 0 : h=min(100*abs(h),h1,hmax)
540 0 : hinit=sign(h,posneg)
541 : return
542 : end function hinit
543 : !
544 :
545 0 : real(dp) function contd5(ii,x,rwork,iwork,ierr)
546 : ! ----------------------------------------------------------
547 : ! this function can be used for continuous output in connection
548 : ! with the output-subroutine for dopri5. it provides an
549 : ! approximation to the ii-th component of the solution at x.
550 : ! ----------------------------------------------------------
551 : integer, intent(in) :: ii ! result is interpolated approximation of y(i) at x=s.
552 : real(dp), intent(in) :: x ! interpolation x value (between xold and x).
553 : real(dp), intent(inout), target :: rwork(*)
554 : integer, intent(inout), target :: iwork(*)
555 : integer, intent(out) :: ierr
556 :
557 0 : real(dp) :: xold, h, theta, theta1
558 : integer :: nd, i, j
559 0 : real(dp), pointer :: con(:)
560 0 : integer, pointer :: icomp(:)
561 :
562 : !dimension con(5*nd),icomp(nd)
563 : !common /condo5/xold,h
564 :
565 0 : nd = iwork(1)
566 0 : icomp => iwork(2:nd+1)
567 0 : xold = rwork(1)
568 0 : h = rwork(2)
569 0 : con => rwork(3:2+5*nd)
570 :
571 : ! ----- compute place of ii-th component
572 0 : i=0
573 0 : do j=1,nd
574 0 : if (icomp(j) == ii) i=j
575 : end do
576 0 : if (i == 0) then
577 0 : contd5 = 0
578 0 : ierr = -1
579 0 : return
580 : end if
581 0 : ierr=0
582 0 : theta=(x-xold)/h
583 0 : theta1=1.d0-theta
584 : contd5=con(i)+theta*(con(nd+i)+theta1*(con(2*nd+i)+theta*
585 0 : & (con(3*nd+i)+theta1*con(4*nd+i))))
586 0 : end function contd5
587 : !
588 :
589 11 : subroutine cdopri(c2,c3,c4,c5,e1,e3,e4,e5,e6,e7,
590 : & a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
591 : & a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,
592 : & d1,d3,d4,d5,d6,d7)
593 : ! ----------------------------------------------------------
594 : ! runge-kutta coefficients of dormand and prince (1980)
595 : ! ----------------------------------------------------------
596 : real(dp) :: c2, c3, c4, c5
597 : real(dp) :: a21
598 : real(dp) :: a31, a32
599 : real(dp) :: a41, a42, a43
600 : real(dp) :: a51, a52, a53, a54
601 : real(dp) :: a61, a62, a63, a64, a65
602 : real(dp) :: a71, a72, a73, a74, a75, a76
603 : real(dp) :: e1, e2, e3, e4, e5, e6, e7
604 : real(dp) :: d1, d3, d4, d5, d6, d7
605 11 : c2=0.2d0
606 11 : c3=0.3d0
607 11 : c4=0.8d0
608 11 : c5=8.d0/9.d0
609 11 : a21=0.2d0
610 11 : a31=3.d0/40.d0
611 11 : a32=9.d0/40.d0
612 11 : a41=44.d0/45.d0
613 11 : a42=-56.d0/15.d0
614 11 : a43=32.d0/9.d0
615 11 : a51=19372.d0/6561.d0
616 11 : a52=-25360.d0/2187.d0
617 11 : a53=64448.d0/6561.d0
618 11 : a54=-212.d0/729.d0
619 11 : a61=9017.d0/3168.d0
620 11 : a62=-355.d0/33.d0
621 11 : a63=46732.d0/5247.d0
622 11 : a64=49.d0/176.d0
623 11 : a65=-5103.d0/18656.d0
624 11 : a71=35.d0/384.d0
625 11 : a73=500.d0/1113.d0
626 11 : a74=125.d0/192.d0
627 11 : a75=-2187.d0/6784.d0
628 11 : a76=11.d0/84.d0
629 11 : e1=71.d0/57600.d0
630 11 : e3=-71.d0/16695.d0
631 11 : e4=71.d0/1920.d0
632 11 : e5=-17253.d0/339200.d0
633 11 : e6=22.d0/525.d0
634 11 : e7=-1.d0/40.d0
635 : ! ---- dense output of shampine (1986)
636 11 : d1=-12715105075.d0/11282082432.d0
637 11 : d3=87487479700.d0/32700410799.d0
638 11 : d4=-10690763975.d0/1880347072.d0
639 11 : d5=701980252875.d0/199316789632.d0
640 11 : d6=-1453857185.d0/822651844.d0
641 11 : d7=69997945.d0/29380423.d0
642 11 : return
643 : end subroutine cdopri
644 :
645 : end module mod_dopri5
|
