EFTCAMB  Reference documentation for version 3.0
02_root_finding.f90
Go to the documentation of this file.
1 !----------------------------------------------------------------------------------------
2 !
3 ! This file is part of EFTCAMB.
4 !
5 ! Copyright (C) 2013-2016 by the EFTCAMB authors
6 !
7 ! The EFTCAMB code is free software;
8 ! You can use it, redistribute it, and/or modify it under the terms
9 ! of the GNU General Public License as published by the Free Software Foundation;
10 ! either version 3 of the License, or (at your option) any later version.
11 ! The full text of the license can be found in the file eftcamb/LICENSE at
12 ! the top level of the EFTCAMB distribution.
13 !
14 !----------------------------------------------------------------------------------------
15 
18 
19 
20 !----------------------------------------------------------------------------------------
22 
24 
26 
27  use precision
28 
29  implicit none
30 
31  private
32 
33  public zbrac, zbrent
34 
35 contains
36 
37  ! ---------------------------------------------------------------------------------------------
40  subroutine zbrac( func, x1, x2, succes, funcZero )
41 
42  implicit none
43 
44  real(dl) func
45  real(dl), intent(inout) :: x1
47  real(dl), intent(inout) :: x2
49  logical , intent(out) :: succes
50  real(dl), intent(in) :: funcZero
51 
52  integer , parameter :: ntry = 1000
53  real(dl), parameter :: FACTOR = 5._dl
54 
55  real(dl) delta, temp, f1,f2
56  integer j
57  external func
58 
59  ! initial check:
60  if ( x1.eq.x2 ) stop 'you have to guess an initial range in zbrac'
61  ! get an ordered interval:
62  if (x2<x1) then
63  temp=x2
64  x2=x1
65  x1=temp
66  end if
67  ! compute the function:
68  f1 = func(x1) -funczero
69  f2 = func(x2) -funczero
70 
71  succes=.true.
72  do j=1, ntry
73  ! check wether the interval is bracketing:
74  if ( f1*f2 .lt. 0._dl ) return
75  ! compute the interval and enlarge it:
76  delta=abs(x2-x1)
77  x1 = x1 -factor*delta
78  x2 = x2 +factor*delta
79  ! recompute the functions:
80  f1 = func(x1) -funczero
81  f2 = func(x2) -funczero
82  end do
83  succes=.false.
84 
85  return
86 
87  end subroutine zbrac
88 
89  ! ---------------------------------------------------------------------------------------------
93  function zbrent(func,x1,x2,tol,funcZero,succes)
94 
95  implicit none
96 
97  real(dl) :: func
98  real(dl) :: x1
99  real(dl) :: x2
100  real(dl) :: tol
102  real(dl) :: funcZero
103  logical :: succes
104 
105  real(dl) :: zbrent
106 
107  external func
108 
109  integer ,parameter :: ITMAX = 2000
110  real(dl),parameter :: EPS = 3.d-15
111  integer iter
112  real(dl) a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
113 
114  succes = .true.
115  a=x1
116  b=x2
117  fa=func(a)-funczero
118  fb=func(b)-funczero
119  if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) then
120  succes=.false.
121  return
122  end if
123  c=b
124  fc=fb
125  do 11 iter=1,itmax
126  if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
127  c=a
128  fc=fa
129  d=b-a
130  e=d
131  endif
132  if(abs(fc).lt.abs(fb)) then
133  a=b
134  b=c
135  c=a
136  fa=fb
137  fb=fc
138  fc=fa
139  endif
140  tol1=2.*eps*abs(b)+0.5*tol
141  xm=.5*(c-b)
142  if(abs(xm).le.tol1 .or. fb.eq.0.)then
143  zbrent=b
144  return
145  endif
146  if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
147  s=fb/fa
148  if(a.eq.c) then
149  p=2.*xm*s
150  q=1.-s
151  else
152  q=fa/fc
153  r=fb/fc
154  p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
155  q=(q-1.)*(r-1.)*(s-1.)
156  endif
157  if(p.gt.0.) q=-q
158  p=abs(p)
159  if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
160  e=d
161  d=p/q
162  else
163  d=xm
164  e=d
165  endif
166  else
167  d=xm
168  e=d
169  endif
170  a=b
171  fa=fb
172  if(abs(d) .gt. tol1) then
173  b=b+d
174  else
175  b=b+sign(tol1,xm)
176  endif
177  fb=func(b)-funczero
178 11 continue
179  succes = .false.
180  zbrent=b
181  return
182 
183  end function zbrent
184 
185 end module eftcamb_rootfind
186 
187 !----------------------------------------------------------------------------------------
real(dl) function, public zbrent(func, x1, x2, tol, funcZero, succes)
Brent root finding algorithm: This is used to solve numerically the equation func=funcZero by means o...
This module contains the definitions of all the EFTCAMB root finding algorithms.
subroutine, public zbrac(func, x1, x2, succes, funcZero)
Bracketing subroutine: This subroutine does a outward search for the smallest intervall containing a ...