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 :
|
