Branch data Line data Source code
1 : : !!****m* ABINIT/m_brentq
2 : : !! NAME
3 : : !! m_brentq
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains Brent's root-finding method.
7 : : !! This was translated from scipy.
8 : : !! Originally written by Charles Harris charles.harris@sdl.usu.edu
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (c) 2001-2002 Enthought, Inc. 2003-2024, SciPy Developers.
12 : : !! All rights reserved.
13 : : !!
14 : : !! Redistribution and use in source and binary forms, with or without
15 : : !! modification, are permitted provided that the following conditions
16 : : !! are met:
17 : : !!
18 : : !! 1. Redistributions of source code must retain the above copyright
19 : : !! notice, this list of conditions and the following disclaimer.
20 : : !!
21 : : !! 2. Redistributions in binary form must reproduce the above
22 : : !! copyright notice, this list of conditions and the following
23 : : !! disclaimer in the documentation and/or other materials provided
24 : : !! with the distribution.
25 : : !!
26 : : !! 3. Neither the name of the copyright holder nor the names of its
27 : : !! contributors may be used to endorse or promote products derived
28 : : !! from this software without specific prior written permission.
29 : : !!
30 : : !! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
31 : : !! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
32 : : !! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
33 : : !! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
34 : : !! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
35 : : !! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
36 : : !! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
37 : : !! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
38 : : !! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
39 : : !! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
40 : : !! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
41 : : !!
42 : : !! This file is distributed under the terms of the
43 : : !! GNU General Public License, see ~abinit/COPYING
44 : : !! or http://www.gnu.org/copyleft/gpl.txt .
45 : : !!
46 : : !! SOURCE
47 : :
48 : : #if defined HAVE_CONFIG_H
49 : : #include "config.h"
50 : : #endif
51 : :
52 : : #include "abi_common.h"
53 : :
54 : : MODULE m_brentq
55 : :
56 : : use defs_basis
57 : : use m_errors
58 : :
59 : : implicit none
60 : :
61 : : private
62 : :
63 : : public :: brentq
64 : :
65 : : CONTAINS !========================================================================================
66 : : !!***
67 : :
68 : : !!****f* m_brentq/brentq
69 : : !! NAME
70 : : !! brentq
71 : : !!
72 : : !! FUNCTION
73 : : !!
74 : : !! Brent's root finding method (taken from scipy).
75 : : !!
76 : : !! INPUTS
77 : : !! f= subroutine corresponding to the function for which the optimization is to be done, of the form subroutine f(x,fx)
78 : : !! xa,xb= interval on which to perform the optimization
79 : : !! xtol,rtol=tolerance criterion, the output root x0 will satisfy the criterion
80 : : !! abs(x-x0) < xtol + rtol * x0 with x the true root, default value on scipy are
81 : : !! 2e-12 and 4*machine_precision respectively
82 : : !! iter= maximum number of iterations, default value on scipy is 100
83 : : !!
84 : : !! OUTPUT
85 : : !! xcur= root
86 : : !! ierr= 1 if the algorithm converged, 0 otherwise
87 : : !!
88 : : !! SOURCE
89 : :
90 : 0 : subroutine brentq(f,xa,xb,xtol,rtol,iter,xcur,ierr)
91 : :
92 : : ! Written by Charles Harris charles.harris@sdl.usu.edu
93 : :
94 : : !Arguments ------------------------------------
95 : : interface
96 : : subroutine f(x,fx)
97 : : use defs_basis
98 : : real(dp), intent(in) :: x
99 : : real(dp), intent(out) :: fx
100 : : end subroutine f
101 : : end interface
102 : : real(dp), intent(in) :: xa,xb,xtol,rtol
103 : : integer, intent(in) :: iter
104 : : real(dp), intent(out) :: xcur
105 : : integer, intent(out) :: ierr
106 : : !Local variables ------------------------------
107 : : integer :: i
108 : : real(dp) :: dblk,delta,dpre,fblk,fcur,fpre,sbis,scur,spre,stry,xblk,xpre
109 : : character(len=200) :: msg
110 : : !************************************************************************
111 : :
112 : 0 : xpre = xa
113 : 0 : xcur = xb
114 : 0 : xblk = zero
115 : 0 : fblk = zero
116 : 0 : spre = zero
117 : 0 : scur = zero
118 : 0 : ierr = 0
119 : :
120 : : ! the tolerance is 2*delta
121 : 0 : call f(xpre,fpre)
122 : 0 : call f(xcur,fcur)
123 : :
124 [ # # ]: 0 : if (fpre == zero) then
125 : 0 : ierr = 1
126 : 0 : xcur = xpre
127 : 0 : return
128 : : end if
129 : :
130 [ # # ]: 0 : if (fcur == zero) then
131 : 0 : ierr = 1
132 : 0 : return
133 : : end if
134 : :
135 [ # # ]: 0 : if (sign(one,fpre) == sign(one,fcur)) then
136 : 0 : msg = 'Sign error in brentq'
137 : 0 : ABI_ERROR(msg)
138 : : end if
139 : :
140 [ # # ]: 0 : do i=1,iter
141 : :
142 [ # # # # : 0 : if ((fpre /= zero) .and. (fcur /= zero) .and. (sign(one,fpre) /= sign(one,fcur))) then
# # ]
143 : 0 : xblk = xpre
144 : 0 : fblk = fpre
145 : 0 : spre = xcur - xpre
146 : 0 : scur = xcur - xpre
147 : : end if
148 : :
149 [ # # ]: 0 : if (abs(fblk) < abs(fcur)) then
150 : 0 : xpre = xcur
151 : 0 : xcur = xblk
152 : 0 : xblk = xpre
153 : 0 : fpre = fcur
154 : 0 : fcur = fblk
155 : 0 : fblk = fpre
156 : : end if
157 : :
158 : 0 : delta = (xtol+rtol*abs(xcur)) * half
159 : 0 : sbis = (xblk - xcur) * half
160 [ # # # # ]: 0 : if ((fcur == zero) .or. (abs(sbis) < delta)) then
161 : 0 : ierr = 1
162 : 0 : return
163 : : end if
164 : :
165 [ # # # # ]: 0 : if ((abs(spre) > delta) .and. (abs(fcur) < abs(fpre))) then
166 [ # # ]: 0 : if (xpre == xblk) then
167 : : ! interpolate
168 : 0 : stry = -fcur * (xcur-xpre) / (fcur-fpre)
169 : : else
170 : : ! extrapolate
171 : 0 : dpre = (fpre-fcur) / (xpre-xcur)
172 : 0 : dblk = (fblk-fcur) / (xblk-xcur)
173 : 0 : stry = -fcur * (fblk*dblk-fpre*dpre) / (dblk*dpre*(fblk-fpre))
174 : : end if
175 : :
176 [ # # ]: 0 : if (two*abs(stry) < min(abs(spre),three*abs(sbis)-delta)) then
177 : : ! good short step
178 : : spre = scur
179 : : scur = stry
180 : : else
181 : : ! bisect
182 : 0 : spre = sbis
183 : 0 : scur = sbis
184 : : end if
185 : : else
186 : : ! bisect
187 : : spre = sbis
188 : : scur = sbis
189 : : end if
190 : :
191 : 0 : xpre = xcur
192 : 0 : fpre = fcur
193 : :
194 [ # # ]: 0 : if (abs(scur) > delta) then
195 : 0 : xcur = xcur + scur
196 : : else
197 [ # # ]: 0 : if (sbis > zero) then
198 : 0 : xcur = xcur + delta
199 : : else
200 : 0 : xcur = xcur - delta
201 : : end if
202 : : end if
203 : :
204 : 0 : call f(xcur,fcur)
205 : :
206 : : end do ! i
207 : :
208 : : end subroutine brentq
209 : : !!***
210 : :
211 : : !----------------------------------------------------------------------
212 : :
213 : : END MODULE m_brentq
214 : : !!***
215 : :
|
