TABLE OF CONTENTS

ABINIT/brent [ Functions ]

[ Top ] [ Functions ]

NAME

 brent

FUNCTION

 minimizes a function along a line

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~ABINIT/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~ABINIT/Infos/contributors .

INPUTS

 dp_dum_vdp: function  to be minimized (return a dp from a vector of dp)
 vdp_dum_vdp: derivative of the function (return a vector of dp from a vector of dp)
 itmax: number of iterations allowed
 tol: tolerance on error. It depend on the precision of the numbers
 (usualy chosen as sqrt(max precision available with your floating point reresentation))
 ax,xx,bx: a bracketing triplet around the minimum to be find

OUTPUT

 xmin: value such that dp_dum_vdp(v(:)+xmin*grad(:)) is minimum
 brent:  dp_dum_vdp(v(:)+xmin*grad(:))

SIDE EFFECTS

 grad(:): direction along which the minimization is performed
 v(:): starting and ending point of the minimization

PARENTS

 linmin

CHILDREN

 dotproduct

SOURCE

 39 #if defined HAVE_CONFIG_H
 40 #include "config.h"
 41 #endif
 42 
 43 #include "abi_common.h"
 44 
 45 
 46 function brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,itmax,v,grad,ax,xx,bx,tol,xmin)
 47 
 48  use defs_basis
 49 
 50 !This section has been created automatically by the script Abilint (TD).
 51 !Do not modify the following lines by hand.
 52 #undef ABI_FUNC
 53 #define ABI_FUNC 'brent'
 54  use interfaces_62_cg_noabirule, except_this_one => brent
 55 !End of the abilint section
 56 
 57  implicit none
 58 
 59 !Arguments ------------------------------------
 60  include "dummy_functions.inc"
 61 !scalars
 62  integer,intent(in) :: itmax,nv1,nv2
 63  real(dp) :: brent
 64  real(dp),intent(in) :: ax,bx,tol,xx
 65  real(dp),intent(out) :: xmin
 66 !arrays
 67  real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2)
 68 
 69 !Local variables-------------------------------
 70 !scalars
 71  integer :: iter
 72  real(dp) :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,vv,w
 73  real(dp) :: x,xm,zeps
 74  logical :: ok1,ok2,ok3,ok4
 75 
 76 !************************************************************************
 77  zeps=epsilon(ax*real(1e-2,dp))
 78  a=min(ax,bx)
 79  b=max(ax,bx)
 80  vv=xx
 81  w=xx
 82  x=xx
 83  e=zero
 84  fx=dp_dum_v2dp(nv1,nv2,x*grad(:,:)+v(:,:))
 85  fv=fx
 86  fw=fx
 87 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of
 88 !v(:,:)=v(:,:)+(grad(:,:)*x)
 89 !but for instance renormilizing the density if brent is used on a density...
 90 !vp(:,:) = v(:,:)
 91 !sub_dum_dp_v2dp_v2dp(x,grad(:,:),vp(:,:)
 92 !dx=dotproduct(v2dp_dum_v2dp(vp(:,:)),grad(:,:))
 93  dx=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,v(:,:)+x*grad(:,:)),grad(:,:))
 94  dv=dx
 95  dw=dx
 96  do iter=1,itmax
 97   xm=half*(a+b)
 98   tol1=tol*abs(x)+zeps
 99   tol2=two*tol1
100   if(abs(x-xm) <= (tol2-half*(b-a))) then
101    exit
102   end if
103   if(abs(e) > tol1) then
104    d1=two*(b-a)
105    d2=d1
106    if(dw /= dx) d1=(w-x)*dx/(dx-dw)
107    if(dv /= dx) d2=(vv-x)*dx/(dx-dv)
108    u1=x+d1
109    u2=x+d2
110    ok1=((a-u1)*(u1-b)>zero).and.(dx*d1<=zero)
111    ok2=((a-u2)*(u2-b)>zero).and.(dx*d2<=zero)
112    olde=e
113    e=d
114    if(ok1.or.ok2) then
115     if(ok1.and.ok2) then
116      d=merge(d1,d2,abs(d1)<abs(d2))
117     else
118      d=merge(d1,d2,ok1)
119     end if
120     if(abs(d)<=abs(half*olde)) then
121      u=x+d
122      if(((u-a)<tol2).or.((b-u)<tol2)) d=sign(tol1,xm-x)
123     else
124      e=merge(a,b,dx>=zero)-x
125      d=half*e
126     end if
127    else
128     e=merge(a,b,dx>=zero)-x
129     d=half*e
130    end if
131   else
132    e=merge(a,b,dx>=zero)-x
133    d=half*e
134   end if
135 
136   if(abs(d) >=tol1)then
137    u=x+d
138    fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
139   else
140    u=x+sign(tol1,d)
141    fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
142    if(fu>fx) then
143     exit
144    end if
145   end if
146   du=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)),grad(:,:))
147   if(fu<=fx)then
148    if(u>=x)then
149     a=x
150    else
151     b=x
152    end if
153    vv=w
154    fv=fw
155    dv=dw
156    w=x
157    fw=fx
158    dw=dx
159    x=u
160    dx=du
161    fx=fu
162   else
163    if(u<x) then
164     a=u
165    else
166     b=u
167    end if
168    ok3=(w==x).or.(fu.le.fw)
169    ok4=(vv==w).or.(vv==x).or.(fu.lt.fv)
170    if(ok3) then
171     vv=w
172     fv=fw
173     dv=dw
174     w=u
175     fw=fu
176     dw=du
177    else if( ok4 ) then
178     vv=u
179     fv=fu
180     dv=du
181    end if
182   end if
183  end do
184  xmin=x
185 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of
186 !v(:,:)=v(:,:)+(grad(:,:)*x)
187 !but for instance renormilizing the density if brent is used on a density...
188  call sub_dum_dp_v2dp_v2dp(nv1,nv2,x,grad(:,:),v(:,:))
189  brent=fx
190 end function brent