https://github.com/cran/nleqslv
Tip revision: 373fd6f8a8ac368501e3ecb1782ca64ca749eda8 authored by Berend Hasselman on 22 April 2009, 16:02:21 UTC
version 1.4
version 1.4
Tip revision: 373fd6f
nwtrup.f
subroutine nwtrup(n,fcnorm,g,sc,scalex,nwtake,stepmx,xtol,dlt,
* mxtake, fpred,retcd,xprev,fpnsav,fprev,xp,fp,
* fpnorm,wrk)
integer n,retcd
double precision fcnorm,stepmx,xtol,dlt,fpred,fpnsav,fpnorm
double precision xp(*),g(*),scalex(*)
double precision sc(*),xprev(*),fprev(*),fp(*),wrk(*)
logical nwtake,mxtake
c-------------------------------------------------------------------------
c
c Decide whether to accept xp=xc+sc as the next iterate
c and updates the trust region dlt
c
c Arguments
c
c In n Integer size of xc()
c In fcnorm Real .5*||f(xc)||**2
c In g Real(*) gradient at xc()
c In sc Real(*) current step
c In scalex Real(*) scaling factors x(*)
c In nwtake Logical true if sc is newton direction
c In stepmx Real maximum step size
c In xtol Real minimum step tolerance
c Inout dlt Real trust region radius
c In mxtake Logical true if max. step taken
c In fpred Real predicted value of .5*||f()||**2
c
c Inout retcd Integer return code
c 0 xp accepted as next iterate;
c dlt trust region for next iteration.
c
c 1 xp unsatisfactory but
c accepted as next iterate because
c xp-xc .lt. smallest allowable
c step length.
c
c 2 f(xp) too large.
c continue current iteration with
c new reduced dlt.
c
c 3 f(xp) sufficiently small, but
c quadratic model predicts f(xp)
c sufficiently well to continue current
c iteration with new doubled dlt.
c
c On first entry, retcd must be 4
c
c Wk xprev Real(*) (internal) work
c Wk fpnsav Real (internal) work
c Wk fprev Real(*) (internal) work
c Inout xp Real(*) new iterate x()
c Inout fp Real(*) new f(xp)
c Inout fpnorm Real new .5*||f(xp)||**2
c Wk wrk Real(*) (internal) work
c
c-------------------------------------------------------------------------
double precision ared,pred,slope,sclen,rln,dltmp
double precision dnrm2,ddot,nudnrm
logical ret3ok
double precision Rone, Rtwo, Rthree, Rfour, Rten
double precision Rhalf, Rpten, Rphund
parameter(Rpten = 0.1d0)
parameter(Rhalf=0.5d0)
parameter(Rone=1.0d0, Rtwo=2.0d0, Rthree=3.0d0, Rfour=4.0d0)
parameter(Rten=10.0d0)
double precision Rp99,Rp4, Rp75
parameter(Rp99=Rone-Rten**(-2), Rp4=Rten**(-4), Rp75=Rthree/Rfour)
double precision alpha
parameter(alpha = Rp4)
mxtake = .false.
c ared measures the actual reduction in the function value
c pred measures the predicted reduction in the function value
ared = fpnorm - fcnorm
pred = fpred - fcnorm
slope = ddot(n,g,1,sc,1)
if(retcd .ne. 3) then
ret3ok = .false.
else
ret3ok = fpnorm .ge. fpnsav .or. ared .gt. alpha * slope
endif
if(retcd .eq. 3 .and. ret3ok) then
c reset xp to xprev and terminate global step
retcd = 0
call dcopy(n,xprev,1,xp,1)
call dcopy(n,fprev,1,fp,1)
fpnorm = fpnsav
dlt = Rhalf*dlt
elseif(ared .gt. alpha * slope) then
c fpnorm too large (decrease not sufficient)
rln = nudnrm(n,sc,xp, scalex)
if(rln .lt. xtol) then
c cannot find satisfactory xp sufficiently distinct from xc
retcd = 1
else
c reduce trust region and continue global step
retcd = 2
call dcopy(n,sc,1,wrk,1)
call vscal(n,wrk,scalex)
sclen = dnrm2(n,wrk,1)
dltmp = -slope*sclen/(Rtwo*(ared-slope))
if(dltmp .lt. Rpten*dlt) then
dlt = Rpten*dlt
else
dlt = min(Rhalf*dlt, dltmp)
endif
endif
elseif(retcd .ne. 2 .and. (abs(pred-ared) .le. Rpten*abs(ared))
* .and. (.not.nwtake) .and. (dlt.le. Rp99*stepmx)) then
c fpnorm sufficiently small
c to attempt a doubling of the trust region and continue global step
call dcopy(n,xp,1,xprev,1)
call dcopy(n,fp,1,fprev,1)
fpnsav = fpnorm
dlt = min(Rtwo*dlt,stepmx)
retcd = 3
else
c fpnorm sufficiently small to accept xp as next iterate.
c Choose new trust region.
retcd = 0
if(dlt .gt. Rp99*stepmx) mxtake = .true.
if(ared .ge. Rpten*pred) then
c Not good enough. Decrease trust region for next iteration
dlt = Rhalf*dlt
elseif( ared .le. Rp75*pred ) then
c Wonderful. Increase trust region for next iteration
dlt = min(Rtwo*dlt,stepmx)
endif
endif
return
end