https://github.com/cran/nleqslv
Tip revision: e881fc8a6d3dd7dc7e2787e219afb5bb914ee360 authored by Berend Hasselman on 04 May 2011, 15:31:55 UTC
version 1.8.4
version 1.8.4
Tip revision: e881fc8
nwnwtn.f
subroutine nwsolv(ldr,xc,n,scalex,maxit,
* jacflg,xtol,ftol,btol,global,xscalm,
* stepmx,dlt,sigma,
* rjac,wrk1,wrk2,wrk3,wrk4,fc,fq,dn,d,qtf,
* rcdwrk,icdwrk,qrwork,qrwsiz,
* epsm,fjac,fvec,outopt,xp,fp,gp,njcnt,nfcnt,
* termcd)
integer ldr,n,termcd,njcnt,nfcnt
integer maxit,jacflg,global,xscalm,qrwsiz
integer outopt(*)
double precision xtol,ftol,btol,stepmx,dlt,sigma,fpnorm,epsm
double precision rjac(ldr,*)
double precision xc(*),fc(*),xp(*),fp(*),dn(*),d(*)
double precision wrk1(*),wrk2(*),wrk3(*),wrk4(*)
double precision qtf(*),gp(*),fq(*)
double precision scalex(*)
double precision rcdwrk(*),qrwork(*)
integer icdwrk(*)
external fjac,fvec
c-----------------------------------------------------------------------
c
c Solve system of nonlinear equations with Newton and global strategy
c
c
c Arguments
c
c In ldr Integer leading dimension of rjac
c In xc Real(*) initial estimate of solution
c In n Integer dimensions of problem
c Inout scalex Real(*) scaling factors x(*)
c In maxit Integer maximum number of allowable iterations
c In jacflg Integer jacobian flag
c 1 if analytic jacobian supplied
c 0 if analytic jacobian not supplied
c In xtol Real tolerance at which successive iterates x()
c are considered close enough to
c terminate algorithm
c In ftol Real tolerance at which function values f()
c are considered close enough to zero
c Inout btol Real x tolerance for backtracking
c In global Integer global strategy to use
c 0 quadratic line search
c 1 geometric line search
c 2 double dogleg
c 3 single dogleg
c In xscalm Integer x scaling method
c 1 from column norms of first jacobian
c increased if needed after first iteration
c 0 scaling user supplied
c In stepmx Real maximum allowable step size
c In dlt Real trust region radius
c In sigma Real reduction factor geometric linesearch
c Inout rjac Real(ldr,*) jacobian (2*n columns)
c Wk wrk1 Real(*) workspace
c Wk wrk2 Real(*) workspace
c Wk wrk3 Real(*) workspace
c Wk wrk4 Real(*) workspace
c Inout fc Real(*) function values f(xc)
c Wk fq Real(*) workspace
c Wk dn Real(*) workspace
c Wk d Real(*) workspace
c Wk qtf Real(*) workspace
c Wk rcdwrk Real(*) workspace
c Wk icdwrk Integer(*) workspace
c In qrwork Real(*) workspace for Lapack QR routines (call nwqmem)
c In qrwsiz Integer size of qrwork
c In epsm Real machine precision
c In fjac Name name of routine to calculate jacobian
c (optional)
c In fvec Name name of routine to calculate f()
c In outopt Integer(*) output options
c Out xp Real(*) final x()
c Out fp Real(*) final f(xp)
c Out gp Real(*) gradient at xp()
c Out njcnt Integer number of jacobian evaluations
c Out nfcnt Integer number of function evaluations
c Out termcd Integer termination code
c
c-----------------------------------------------------------------------
integer iter,gcnt,retcd,ierr
double precision dum(2),fcnorm,rcond
logical mxtake
integer priter
integer idamax
double precision Rzero, Rone
parameter(Rzero=0.0d0, Rone=1.0d0)
c initialization
retcd = 0
iter = 0
njcnt = 0
nfcnt = 0
ierr = 0
dum(1) = 0
if( outopt(1) .eq. 1 ) then
priter = 1
else
priter = -1
endif
c evaluate function
call vscal(n,xc,scalex)
call nwfvec(xc,n,scalex,fvec,fc,fcnorm,wrk1)
c evaluate analytic or finite difference jacobian and check analytic
c jacobian, if requested
if(jacflg .eq. 1) then
if( outopt(2) .eq. 1 ) then
call nwfjac(xc,scalex,fc,fq,n,epsm,jacflg,fvec,fjac,rjac,
* ldr,wrk1)
call chkjac(rjac,ldr,xc,fc,n,epsm,scalex,
* fq,wrk1,fvec,termcd)
if(termcd .lt. 0) then
c copy initial values
call dcopy(n,xc,1,xp,1)
call dcopy(n,fc,1,fp,1)
fpnorm = fcnorm
return
endif
endif
endif
c check stopping criteria for input xc
call nwtcvg(xc,fc,xc,xtol,retcd,ftol,iter,maxit,n,ierr, termcd)
if(termcd .gt. 0) then
call dcopy(n,xc,1,xp,1)
call dcopy(n,fc,1,fp,1)
fpnorm = fcnorm
return
endif
if( priter .gt. 0 ) then
dum(1) = fcnorm
dum(2) = abs(fc(idamax(n,fc,1)))
if( global .le. 1 ) then
call nwlsot(iter,-1,dum)
elseif( global .eq. 2 ) then
call nwdgot(iter,-1,dum)
elseif( global .eq. 3 ) then
call nwpwot(iter,-1,dum)
endif
endif
do while( termcd .eq. 0 )
iter = iter + 1
c - evaluate the jacobian at the current iterate xc
c - evaluate the gradient at the current iterate xc
call nwfjac(xc,scalex,fc,fq,n,epsm,jacflg,fvec,fjac,rjac,
* ldr,wrk1)
njcnt = njcnt + 1
c - if requested calculate x scale from jacobian column norms a la Minpack
if( xscalm .eq. 1 ) then
call vunsc(n,xc,scalex)
call nwcpsx(n,rjac,ldr,scalex,epsm,iter)
call vscal(n,xc,scalex)
endif
call nwscjac(n,rjac,ldr,scalex)
c gp = trans(Rjac) * fc
call dgemv('T',n,n,Rone,rjac,ldr,fc,1,Rzero,gp,1)
c - get newton step
call dcopy(n,fc,1,fq,1)
call nwndir(rjac,ldr,rjac(1,n+1),fq,n,epsm,jacflg,
* wrk1,wrk2,wrk3,wrk4,dn,qtf,ierr,rcond,
* rcdwrk,icdwrk,qrwork,qrwsiz)
call nwsnot(0,ierr,rcond)
c - choose the next iterate xp by a global strategy
if( ierr .gt. 0 ) then
c jacobian singular or too ill-conditioned
call dcopy(n,xc,1,xp,1)
call dcopy(n,fc,1,fp,1)
fpnorm = fcnorm
gcnt = 0
if( priter .gt. 0 ) then
call nwjerr(iter)
endif
elseif(global .eq. 0) then
call nwqlsh(n,xc,fcnorm,dn,gp,stepmx,btol,scalex,
* fvec,xp,fp,fpnorm,wrk1,mxtake,retcd,gcnt,
* priter,iter)
elseif(global .eq. 1) then
call nwglsh(n,xc,fcnorm,dn,gp,sigma,stepmx,btol,scalex,
* fvec,xp,fp,fpnorm,wrk1,mxtake,retcd,gcnt,
* priter,iter)
elseif(global .eq. 2) then
call nwddlg(n,rjac(1,n+1),ldr,dn,gp,xc,fcnorm,stepmx,
* btol,mxtake,dlt,qtf,scalex,
* fvec,d,fq,wrk1,wrk2,wrk3,wrk4,
* xp,fp,fpnorm,retcd,gcnt,priter,iter)
elseif(global .eq. 3) then
call nwpdlg(n,rjac(1,n+1),ldr,dn,gp,xc,fcnorm,stepmx,
* btol,mxtake,dlt,qtf,scalex,
* fvec,d,fq,wrk1,wrk2,wrk3,wrk4,
* xp,fp,fpnorm,retcd,gcnt,priter,iter)
endif
nfcnt = nfcnt + gcnt
c - check stopping criteria for the new iterate xp
call nwtcvg(xp,fp,xc,xtol,retcd,ftol,iter,maxit,n,ierr,termcd)
if(termcd .eq. 0) then
c update xc, fc, and fcnorm
call dcopy(n,xp,1,xc,1)
call dcopy(n,fp,1,fc,1)
fcnorm = fpnorm
endif
enddo
call vunsc(n,xp,scalex)
return
end
c-----------------------------------------------------------------------
subroutine nwndir(rjac,ldr,r,fn,n,epsm,jacflg,
* qraux,y,w,wa,dn,qtf,ierr,rcond,
* rcdwrk,icdwrk,qrwork,qrwsiz)
integer ldr,n,ierr,jacflg,qrwsiz
double precision epsm,rjac(ldr,*),r(ldr,*),qraux(*),fn(*)
double precision wa(*),dn(*),y(*),w(*),qtf(*)
double precision rcdwrk(*),qrwork(*)
integer icdwrk(*)
double precision rcond
c-----------------------------------------------------------------------
c
c Calculate the newton direction
c
c Arguments
c
c Inout rjac Real(ldr,*) jacobian matrix at current iterate
c overwritten with QR decomposition
c In ldr Integer leading dimension of rjac
c Out r Real(ldr,*) upper triangular R from QR decomposition
c In fn Real(*) function values at current iterate
c In n Integer dimension of problem
c In epsm Real machine precision
c In jacflg Integer jacobian flag
c 1 for analytic
c 0 for numeric
c used for condition estimate
c Inout qraux Real(*) QR info from liqrdc
c Wk y Real(*) workspace
c Wk w Real(*) workspace
c Wk wa Real(*) workspace
c Out dn Real(*) Newton direction
c Out qtf Real(*) trans(Q)*f()
c Out ierr Integer 0 indicating Jacobian not ill-conditioned or singular
c 1 indicating Jacobian ill-conditioned
c 2 indicating Jacobian completely singular
c Out rcond Real inverse condition matrix
c Wk rcdwrk Real(*) workspace
c Wk icdwrk Integer(*) workspace
c In qrwork Real(*) workspace for Lapack QR routines (call nwqmem)
c In qrwsiz Integer size of qrwork
c
c-----------------------------------------------------------------------
integer j,info
double precision Rzero, Rone
parameter(Rzero=0.0d0, Rone=1.0d0)
c perform a QR factorization of rjac (simple Lapack routine)
c check for singularity or ill conditioning
c form qtf = trans(Q) * fn
c copy upper triangular part of QR to R
call liqrfa(rjac,ldr,n,qraux,qrwork,qrwsiz,ierr)
c check for singularity or ill conditioning
call cndjac(n,rjac,ldr,epsm,rcond,y,rcdwrk,icdwrk,ierr)
if( ierr .ne. 0 ) then
return
endif
c compute qtf = trans(Q)*fn
call dcopy(n,fn,1,qtf,1)
call liqrqt(rjac, ldr, n, qraux, qtf, qrwork, qrwsiz, info)
c copy the upper triangular part of a QR decomposition
c contained in Rjac into R.
do j=1,n
call dcopy(j,rjac(1,j),1,r(1,j),1)
enddo
c solve rjac*dn = -fn
c ==> R*dn = - qtf
call dcopy(n,qtf,1,dn,1)
call dtrsv('U','N','N',n,rjac,ldr,dn,1)
call dscal(n, -Rone, dn, 1)
return
end