istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/equilibrium: rewrite points_tgo, points_ox

Browse Source 
This change adds a bit of documentation and simplifies the two
(internal) subroutines used to find the horizontal tangent points
and the magnetic O/X point.

Using a closure we can avoid explicitly passing parameters (psi0) to
hybrj1. Previously this required a custom `hybrj1mv` subroutine in
fitpack with an identical interface, except for our extra parameter.
This commit is contained in:
Michele Guerini Rocco 2024-08-27 18:19:28 +02:00
parent c44176a505
commit 15a1f866b4
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
3 changed files with 108 additions and 703 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

225
src/equilibrium.f90
View File

@ -69,12 +69,12 @@ module equilibrium
public unset_equil_spline ! Deinitialising internal state
! Members exposed for magsurf_data
public kspl, psi_spline, points_tgo, model
public kspl, psi_spline, find_htg_point
public model, btaxis, btrcen
! Members exposed to gray_core and more
public psia, phitedge
public btaxis, rmaxis, zmaxis, sgnbphi
public btrcen, rcen
public rmaxis, zmaxis, sgnbphi
public rmnm, rmxm, zmnm, zmxm
public zbinf, zbsup
@ -588,7 +588,7 @@ contains
! Search for exact location of the magnetic axis
rax0=data%rax
zax0=data%zax
call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)
call find_ox_point(rax0,zax0,rmaxis,zmaxis,psinoptmp)
write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
@ -599,7 +599,7 @@ contains
select case (ixploc)
case (X_AT_BOTTOM)
call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
call find_ox_point(rbinf,zbinf,r1,z1,psinxptmp)
if(psinxptmp/=-1.0_wp_) then
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
'r', r1, 'z', z1, 'ψ', psinxptmp
@ -608,14 +608,14 @@ contains
zbinf=z1
psinop=psinoptmp
psiant=psinxptmp-psinop
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one)
zbsup=z1
else
ixploc=0
end if
case (X_AT_TOP)
call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info)
call find_ox_point(rbsup,zbsup,r1,z1,psinxptmp)
if(psinxptmp.ne.-1.0_wp_) then
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
'r', r1, 'z', z1, 'ψ', psinxptmp
@ -624,7 +624,7 @@ contains
zbsup=z1
psinop=psinoptmp
psiant=psinxptmp-psinop
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one)
zbinf=z1
else
ixploc=0
@ -634,11 +634,11 @@ contains
psinop=psinoptmp
psiant=one-psinop
! Find upper horizontal tangent point
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one)
zbsup=z1
rbsup=r1
! Find lower horizontal tangent point
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one)
zbinf=z1
rbinf=r1
write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
@ -1276,142 +1276,125 @@ contains
end function tor_curr
subroutine points_ox(rz,zz,rf,zf,psinvf,info)
! Finds the location of the O,X points
subroutine find_ox_point(R0, z0, R1, z1, psi1)
! Given the point (R₀,z₀) as an initial guess, finds
! the exact location (R₁,z₁) where ∇ψ(R₁,z₁) = 0.
! It also returns ψ₁=ψ(R₁,z₁).
!
! Note: this is used to find the magnetic X and O point,
! because both are stationary points for ψ(R,z).
use const_and_precisions, only : comp_eps
use minpack, only : hybrj1
use logger, only : log_error, log_debug
! local constants
integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
! arguments
real(wp_), intent(in) :: rz,zz
real(wp_), intent(out) :: rf,zf,psinvf
integer, intent(out) :: info
! subroutine arguments
real(wp_), intent(in) :: R0, z0
real(wp_), intent(out) :: R1, z1, psi1
! local variables
real(wp_) :: tol
real(wp_), dimension(n) :: xvec,fvec
real(wp_), dimension(lwa) :: wa
real(wp_), dimension(ldfjac,n) :: fjac
integer :: info
real(wp_) :: sol(2), f(2), df(2,2), wa(15)
character(256) :: msg
xvec(1)=rz
xvec(2)=zz
tol = sqrt(comp_eps)
call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info /= 1) then
write (msg, '("O,X coordinates:",2(x,", ",g0.3))') xvec
call log_debug(msg, mod='equilibrium', proc='points_ox')
sol = [R0, z0] ! first guess
call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, &
tol=sqrt(comp_eps), info=info, wa=wa, lwa=15)
if (info /= 1) then
write (msg, '("guess:", 2(x,", ",g0.3))') R0, z0
call log_debug(msg, mod='equilibrium', proc='find_ox_point')
write (msg, '("solution:", 2(x,", ",g0.3))') sol
call log_debug(msg, mod='equilibrium', proc='find_ox_point')
write (msg, '("hybrj1 failed with error ",g0)') info
call log_error(msg, mod='equilibrium', proc='points_ox')
call log_error(msg, mod='equilibrium', proc='find_ox_point')
end if
rf=xvec(1)
zf=xvec(2)
call pol_flux(rf, zf, psinvf)
end subroutine points_ox
R1 = sol(1) ! solution
z1 = sol(2)
call pol_flux(R1, z1, psi1)
contains
subroutine equation(n, x, f, df, ldf, flag)
! The equation to solve: f(R,z) = ∇ψ(R,z) = 0
! subroutine arguments
integer, intent(in) :: n, flag, ldf
real(wp_), intent(in) :: x(n)
real(wp_), intent(inout) :: f(n), df(ldf,n)
if (flag == 1) then
! return f(R,z) = ∇ψ(R,z)
call pol_flux(R=x(1), z=x(2), dpsidr=f(1), dpsidz=f(2))
else
! return ∇f(R,z) = ∇∇ψ(R,z)
call pol_flux(R=x(1), z=x(2), ddpsidrr=df(1,1), &
ddpsidzz=df(2,2), ddpsidrz=df(1,2))
df(2,1) = df(1,2)
end if
end subroutine equation
end subroutine find_ox_point
subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
use logger, only : log_error
! subroutine arguments
integer, intent(in) :: n,iflag,ldfjac
real(wp_), dimension(n), intent(in) :: x
real(wp_), dimension(n), intent(inout) :: fvec
real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
! local variables
real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
character(64) :: msg
select case(iflag)
case(1)
call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz)
fvec(1) = dpsidr
fvec(2) = dpsidz
case(2)
call pol_flux(x(1), x(2), ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz, &
ddpsidrz=ddpsidrz)
fjac(1,1) = ddpsidrr
fjac(1,2) = ddpsidrz
fjac(2,1) = ddpsidrz
fjac(2,2) = ddpsidzz
case default
write (msg, '("invalid iflag: ",g0)')
call log_error(msg, mod='equilibrium', proc='fcnox')
end select
end subroutine fcnox
subroutine points_tgo(rz,zz,rf,zf,psin0,info)
subroutine find_htg_point(R0, z0, R1, z1, psi0)
! Given the point (R₀,z₀) as an initial guess, finds
! the exact location (R₁,z₁) where:
! { ψ(R₁,z₁) = ψ₀
! { ∂ψ/∂R(R₁,z₁) = 0 .
!
! Note: this is used to find the horizontal tangent
! point of the contour ψ(R,z)=ψ₀.
use const_and_precisions, only : comp_eps
use minpack, only : hybrj1mv
use minpack, only : hybrj1
use logger, only : log_error, log_debug
! local constants
integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
! subroutine arguments
real(wp_), intent(in) :: R0, z0, psi0
real(wp_), intent(out) :: R1, z1
! arguments
real(wp_), intent(in) :: rz,zz,psin0
real(wp_), intent(out) :: rf,zf
integer, intent(out) :: info
! local variables
integer :: info
real(wp_) :: sol(2), f(2), df(2,2), wa(15)
character(256) :: msg
! local variables
real(wp_) :: tol
real(wp_), dimension(n) :: xvec,fvec,f0
real(wp_), dimension(lwa) :: wa
real(wp_), dimension(ldfjac,n) :: fjac
sol = [R0, z0] ! first guess
xvec(1)=rz
xvec(2)=zz
f0(1)=psin0
f0(2)=0.0_wp_
tol = sqrt(comp_eps)
call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info /= 1) then
write (msg, '("R,z coordinates:",5(x,g0.3))') xvec, rz, zz, psin0
call log_debug(msg, mod='equilibrium', proc='points_tgo')
write (msg, '("hybrj1mv failed with error ",g0)') info
call log_error(msg, mod='equilibrium', proc='points_tgo')
call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, &
tol=sqrt(comp_eps), info=info, wa=wa, lwa=15)
if (info /= 1) then
write (msg, '("guess:", 2(x,", ",g0.3))') R0, z0
call log_debug(msg, mod='equilibrium', proc='find_htg_point')
write (msg, '("solution:", 2(x,", ",g0.3))') sol
call log_debug(msg, mod='equilibrium', proc='find_htg_point')
write (msg, '("hybrj1 failed with error ",g0)') info
call log_error(msg, mod='equilibrium', proc='find_htg_point')
end if
rf=xvec(1)
zf=xvec(2)
end
R1 = sol(1) ! solution
z1 = sol(2)
contains
subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag)
use logger, only : log_error
subroutine equation(n, x, f, df, ldf, flag)
! The equation to solve: f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] = 0
! subroutine arguments
integer, intent(in) :: n,ldfjac,iflag
real(wp_), dimension(n), intent(in) :: x,f0
real(wp_), dimension(n), intent(inout) :: fvec
real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
! subroutine arguments
integer, intent(in) :: n, ldf, flag
real(wp_), intent(in) :: x(n)
real(wp_), intent(inout) :: f(n), df(ldf, n)
! local variables
real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz
character(64) :: msg
if (flag == 1) then
! return f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R]
call pol_flux(R=x(1), z=x(2), psi_n=f(1), dpsidr=f(2))
f(1) = f(1) - psi0
else
! return ∇f(R,z) = [[∂ψ/∂R, ∂ψ/∂z], [∂²ψ/∂R², ∂²ψ/∂R∂z]]
call pol_flux(R=x(1), z=x(2), dpsidr=df(1,1), dpsidz=df(1,2), &
ddpsidrr=df(2,1), ddpsidrz=df(2,2))
end if
end subroutine equation
select case(iflag)
case(1)
call pol_flux(x(1), x(2), psinv, dpsidr)
fvec(1) = psinv-f0(1)
fvec(2) = dpsidr-f0(2)
case(2)
call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz, &
ddpsidrr=ddpsidrr, ddpsidrz=ddpsidrz)
fjac(1,1) = dpsidr
fjac(1,2) = dpsidz
fjac(2,1) = ddpsidrr
fjac(2,2) = ddpsidrz
case default
write (msg, '("invalid iflag: ",g0)')
call log_error(msg, mod='equilibrium', proc='fcntgo')
end select
end subroutine fcntgo
end subroutine find_htg_point
subroutine unset_equil_spline

8
src/magsurf_data.f90
View File

@ -452,7 +452,7 @@ contains
use logger, only : log_warning
use dierckx, only : profil, sproota
use equilibrium, only : model, frhotor, psi_spline, &
kspl, points_tgo
kspl, find_htg_point
! local constants
integer, parameter :: mest=4
@ -464,7 +464,7 @@ contains
real(wp_), intent(inout) :: rup, zup, rlw, zlw
! local variables
integer :: npoints,np,info,ic,ier,iopt,m
integer :: npoints,np,ic,ier,iopt,m
real(wp_) :: ra,rb,za,zb,th,zc
real(wp_), dimension(mest) :: zeroc
real(wp_), dimension(psi_spline%nknots_x) :: czc
@ -488,8 +488,8 @@ contains
rb=rlw
za=zup
zb=zlw
call points_tgo(ra,za,rup,zup,h,info)
call points_tgo(rb,zb,rlw,zlw,h,info)
call find_htg_point(ra,za,rup,zup,h)
call find_htg_point(rb,zb,rlw,zlw,h)
rcn(1)=rlw
zcn(1)=zlw

578
src/vendor/minpack.f90 vendored
View File

@ -637,584 +637,6 @@ contains
end subroutine hybrj
subroutine hybrj1mv(fcn,n,x,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
use const_and_precisions, only : zero, one
! arguments
integer, intent(in) :: n, ldfjac, lwa
integer, intent(out) :: info
real(wp_), intent(in) :: tol,f0(n)
real(wp_), intent(out) :: wa(lwa)
real(wp_), intent(inout) :: fvec(n), fjac(ldfjac,n), x(n)
! **********
!
! subroutine hybrj1mv
!
! the purpose of hybrj1mv is to find a zero of a system of
! n nonlinear functions in n variables by a modification
! of the powell hybrid method. this is done by using the
! more general nonlinear equation solver hybrjmv. the user
! must provide a subroutine which calculates the functions
! and the jacobian.
!
! the subroutine statement is
!
! subroutine hybrj1mv(fcn,n,x,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
!
! where
!
! fcn is the name of the user-supplied subroutine which
! calculates the functions and the jacobian. fcn must
! be declared in an external statement in the user
! calling program, and should be written as follows.
!
! subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
! integer n,ldfjac,iflag
! real(8) x(n),fvec(n),fjac(ldfjac,n)
! ----------
! if iflag = 1 calculate the functions at x and
! return this vector in fvec. do not alter fjac.
! if iflag = 2 calculate the jacobian at x and
! return this matrix in fjac. do not alter fvec.
! ---------
! return
! end
!
! the value of iflag should not be changed by fcn unless
! the user wants to terminate execution of hybrj1mv.
! in this case set iflag to a negative integer.
!
! n is a positive integer input variable set to the number
! of functions and variables.
!
! x is an array of length n. on input x must contain
! an initial estimate of the solution vector. on output x
! contains the final estimate of the solution vector.
!
! fvec is an output array of length n which contains
! the functions evaluated at the output x.
!
! fjac is an output n by n array which contains the
! orthogonal matrix q produced by the qr factorization
! of the final approximate jacobian.
!
! ldfjac is a positive integer input variable not less than n
! which specifies the leading dimension of the array fjac.
!
! tol is a nonnegative input variable. termination occurs
! when the algorithm estimates that the relative error
! between x and the solution is at most tol.
!
! info is an integer output variable. if the user has
! terminated execution, info is set to the (negative)
! value of iflag. see description of fcn. otherwise,
! info is set as follows.
!
! info = 0 improper input parameters.
!
! info = 1 algorithm estimates that the relative error
! between x and the solution is at most tol.
!
! info = 2 number of calls to fcn with iflag = 1 has
! reached 100*(n+1).
!
! info = 3 tol is too small. no further improvement in
! the approximate solution x is possible.
!
! info = 4 iteration is not making good progress.
!
! wa is a work array of length lwa.
!
! lwa is a positive integer input variable not less than
! (n*(n+13))/2.
!
! subprograms called
!
! user-supplied ...... fcn
!
! minpack-supplied ... hybrjmv
!
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
!
! **********
! local variables
integer :: j, lr, maxfev, mode, nfev, njev, nprint
real(wp_) :: xtol
! parameters
real(wp_), parameter :: factor=1.0e2_wp_
interface
subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
use const_and_precisions, only : wp_
integer, intent(in) :: n,ldfjac,iflag
real(wp_), intent(in) :: x(n),f0(n)
real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n)
end subroutine fcn
end interface
info = 0
!
! check the input parameters for errors.
!
if (n <= 0 .or. ldfjac < n .or. tol < zero &
.or. lwa < (n*(n + 13))/2) return
!
! call hybrjmv.
!
maxfev = 100*(n + 1)
xtol = tol
mode = 2
do j = 1, n
wa(j) = one
end do
nprint = 0
lr = (n*(n + 1))/2
call hybrjmv(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,wa(1),mode, &
factor,nprint,info,nfev,njev,wa(6*n+1),lr,wa(n+1), &
wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
if (info == 5) info = 4
end subroutine hybrj1mv
subroutine hybrjmv(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, &
factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, &
wa3,wa4)
use const_and_precisions, only : zero, one, epsmch=>comp_eps
! arguments
integer, intent(in) :: n, ldfjac, maxfev, mode, nprint, lr
integer, intent(out) :: info, nfev, njev
real(wp_), intent(in) :: xtol, factor, f0(n)
real(wp_), intent(out) :: fvec(n), fjac(ldfjac,n), r(lr), qtf(n), &
wa1(n), wa2(n), wa3(n), wa4(n)
real(wp_), intent(inout) :: x(n), diag(n)
! **********
!
! subroutine hybrj
!
! the purpose of hybrj is to find a zero of a system of
! n nonlinear functions in n variables by a modification
! of the powell hybrid method. the user must provide a
! subroutine which calculates the functions and the jacobian.
!
! the subroutine statement is
!
! subroutine hybrj(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,diag,
! mode,factor,nprint,info,nfev,njev,r,lr,qtf,
! wa1,wa2,wa3,wa4)
!
! where
!
! fcn is the name of the user-supplied subroutine which
! calculates the functions and the jacobian. fcn must
! be declared in an external statement in the user
! calling program, and should be written as follows.
!
! subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
! integer n,ldfjac,iflag
! real(8) x(n),f0(n),fvec(n),fjac(ldfjac,n)
! ----------
! if iflag = 1 calculate the functions at x and
! return this vector in fvec. do not alter fjac.
! if iflag = 2 calculate the jacobian at x and
! return this matrix in fjac. do not alter fvec.
! ---------
! return
! end
!
! the value of iflag should not be changed by fcn unless
! the user wants to terminate execution of hybrj.
! in this case set iflag to a negative integer.
!
! n is a positive integer input variable set to the number
! of functions and variables.
!
! x is an array of length n. on input x must contain
! an initial estimate of the solution vector. on output x
! contains the final estimate of the solution vector.
!
! fvec is an output array of length n which contains
! the functions evaluated at the output x.
!
! fjac is an output n by n array which contains the
! orthogonal matrix q produced by the qr factorization
! of the final approximate jacobian.
!
! ldfjac is a positive integer input variable not less than n
! which specifies the leading dimension of the array fjac.
!
! xtol is a nonnegative input variable. termination
! occurs when the relative error between two consecutive
! iterates is at most xtol.
!
! maxfev is a positive integer input variable. termination
! occurs when the number of calls to fcn with iflag = 1
! has reached maxfev.
!
! diag is an array of length n. if mode = 1 (see
! below), diag is internally set. if mode = 2, diag
! must contain positive entries that serve as
! multiplicative scale factors for the variables.
!
! mode is an integer input variable. if mode = 1, the
! variables will be scaled internally. if mode = 2,
! the scaling is specified by the input diag. other
! values of mode are equivalent to mode = 1.
!
! factor is a positive input variable used in determining the
! initial step bound. this bound is set to the product of
! factor and the euclidean norm of diag*x if nonzero, or else
! to factor itself. in most cases factor should lie in the
! interval (.1,100.). 100. is a generally recommended value.
!
! nprint is an integer input variable that enables controlled
! printing of iterates if it is positive. in this case,
! fcn is called with iflag = 0 at the beginning of the first
! iteration and every nprint iterations thereafter and
! immediately prior to return, with x and fvec available
! for printing. fvec and fjac should not be altered.
! if nprint is not positive, no special calls of fcn
! with iflag = 0 are made.
!
! info is an integer output variable. if the user has
! terminated execution, info is set to the (negative)
! value of iflag. see description of fcn. otherwise,
! info is set as follows.
!
! info = 0 improper input parameters.
!
! info = 1 relative error between two consecutive iterates
! is at most xtol.
!
! info = 2 number of calls to fcn with iflag = 1 has
! reached maxfev.
!
! info = 3 xtol is too small. no further improvement in
! the approximate solution x is possible.
!
! info = 4 iteration is not making good progress, as
! measured by the improvement from the last
! five jacobian evaluations.
!
! info = 5 iteration is not making good progress, as
! measured by the improvement from the last
! ten iterations.
!
! nfev is an integer output variable set to the number of
! calls to fcn with iflag = 1.
!
! njev is an integer output variable set to the number of
! calls to fcn with iflag = 2.
!
! r is an output array of length lr which contains the
! upper triangular matrix produced by the qr factorization
! of the final approximate jacobian, stored rowwise.
!
! lr is a positive integer input variable not less than
! (n*(n+1))/2.
!
! qtf is an output array of length n which contains
! the vector (q transpose)*fvec.
!
! wa1, wa2, wa3, and wa4 are work arrays of length n.
!
! subprograms called
!
! user-supplied ...... fcn
!
! minpack-supplied ... dogleg,enorm,
! qform,qrfac,r1mpyq,r1updt
!
! fortran-supplied ... abs,dmax1,dmin1,mod
!
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
!
! **********
! local variables
integer :: i, iflag, iter, j, jm1, l, ncfail, ncsuc, nslow1, nslow2
integer, dimension(1) :: iwa
logical :: jeval, sing
real(wp_) :: actred, delta, fnorm, fnorm1, pnorm, prered, &
ratio, summ, temp, xnorm
! parameters
real(wp_), parameter :: p1 = 1.0e-1_wp_, p5 = 5.0e-1_wp_, &
p001 = 1.0e-3_wp_, p0001 = 1.0e-4_wp_
interface
subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
use const_and_precisions, only : wp_
integer, intent(in) :: n,ldfjac,iflag
real(wp_), intent(in) :: x(n),f0(n)
real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n)
end subroutine fcn
end interface
!
info = 0
iflag = 0
nfev = 0
njev = 0
!
! check the input parameters for errors.
!
if (n <= 0 .or. ldfjac < n .or. xtol < zero &
.or. maxfev <= 0 .or. factor <= zero &
.or. lr < (n*(n + 1))/2) go to 300
if (mode == 2) then
do j = 1, n
if (diag(j) <= zero) go to 300
end do
end if
!
! evaluate the function at the starting point
! and calculate its norm.
!
iflag = 1
call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
nfev = 1
if (iflag < 0) go to 300
fnorm = enorm(n,fvec)
!
! initialize iteration counter and monitors.
!
iter = 1
ncsuc = 0
ncfail = 0
nslow1 = 0
nslow2 = 0
!
! beginning of the outer loop.
!
do
jeval = .true.
!
! calculate the jacobian matrix.
!
iflag = 2
call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
njev = njev + 1
if (iflag < 0) go to 300
!
! compute the qr factorization of the jacobian.
!
call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
!
! on the first iteration and if mode is 1, scale according
! to the norms of the columns of the initial jacobian.
!
if (iter == 1) then
if (mode /= 2) then
do j = 1, n
diag(j) = wa2(j)
if (wa2(j) == zero) diag(j) = one
end do
end if
!
! on the first iteration, calculate the norm of the scaled x
! and initialize the step bound delta.
!
do j = 1, n
wa3(j) = diag(j)*x(j)
end do
xnorm = enorm(n,wa3)
delta = factor*xnorm
if (delta == zero) delta = factor
end if
!
! form (q transpose)*fvec and store in qtf.
!
do i = 1, n
qtf(i) = fvec(i)
end do
do j = 1, n
if (fjac(j,j) /= zero) then
summ = zero
do i = j, n
summ = summ + fjac(i,j)*qtf(i)
end do
temp = -summ/fjac(j,j)
do i = j, n
qtf(i) = qtf(i) + fjac(i,j)*temp
end do
end if
end do
!
! copy the triangular factor of the qr factorization into r.
!
sing = .false.
do j = 1, n
l = j
jm1 = j - 1
do i = 1, jm1
r(l) = fjac(i,j)
l = l + n - i
end do
r(l) = wa1(j)
if (wa1(j) == zero) sing = .true.
end do
!
! accumulate the orthogonal factor in fjac.
!
call qform(n,n,fjac,ldfjac,wa1)
!
! rescale if necessary.
!
if (mode /= 2) then
do j = 1, n
diag(j) = dmax1(diag(j),wa2(j))
end do
end if
!
! beginning of the inner loop.
!
do
!
! if requested, call fcn to enable printing of iterates.
!
if (nprint > 0) then
iflag = 0
if (mod(iter-1,nprint) == 0) call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
if (iflag < 0) go to 300
end if
!
! determine the direction p.
!
call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
!
! store the direction p and x + p. calculate the norm of p.
!
do j = 1, n
wa1(j) = -wa1(j)
wa2(j) = x(j) + wa1(j)
wa3(j) = diag(j)*wa1(j)
end do
pnorm = enorm(n,wa3)
!
! on the first iteration, adjust the initial step bound.
!
if (iter == 1) delta = dmin1(delta,pnorm)
!
! evaluate the function at x + p and calculate its norm.
!
iflag = 1
call fcn(n,wa2,f0,wa4,fjac,ldfjac,iflag)
nfev = nfev + 1
if (iflag < 0) go to 300
fnorm1 = enorm(n,wa4)
!
! compute the scaled actual reduction.
!
actred = -one
if (fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
!
! compute the scaled predicted reduction.
!
l = 1
do i = 1, n
summ = zero
do j = i, n
summ = summ + r(l)*wa1(j)
l = l + 1
end do
wa3(i) = qtf(i) + summ
end do
temp = enorm(n,wa3)
prered = zero
if (temp < fnorm) prered = one - (temp/fnorm)**2
!
! compute the ratio of the actual to the predicted
! reduction.
!
ratio = zero
if (prered > zero) ratio = actred/prered
!
! update the step bound.
!
if (ratio < p1) then
ncsuc = 0
ncfail = ncfail + 1
delta = p5*delta
else
ncfail = 0
ncsuc = ncsuc + 1
if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5)
if (abs(ratio-one) <= p1) delta = pnorm/p5
end if
!
! test for successful iteration.
!
if (ratio >= p0001) then
!
! successful iteration. update x, fvec, and their norms.
!
do j = 1, n
x(j) = wa2(j)
wa2(j) = diag(j)*x(j)
fvec(j) = wa4(j)
end do
xnorm = enorm(n,wa2)
fnorm = fnorm1
iter = iter + 1
end if
!
! determine the progress of the iteration.
!
nslow1 = nslow1 + 1
if (actred >= p001) nslow1 = 0
if (jeval) nslow2 = nslow2 + 1
if (actred >= p1) nslow2 = 0
!
! test for convergence.
!
if (delta <= xtol*xnorm .or. fnorm == zero) info = 1
if (info /= 0) go to 300
!
! tests for termination and stringent tolerances.
!
if (nfev >= maxfev) info = 2
if (p1*dmax1(p1*delta,pnorm) <= epsmch*xnorm) info = 3
if (nslow2 == 5) info = 4
if (nslow1 == 10) info = 5
if (info /= 0) go to 300
!
! criterion for recalculating jacobian.
!
if (ncfail == 2) exit
!
! calculate the rank one modification to the jacobian
! and update qtf if necessary.
!
do j = 1, n
summ = zero
do i = 1, n
summ = summ + fjac(i,j)*wa4(i)
end do
wa2(j) = (summ - wa3(j))/pnorm
wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
if (ratio >= p0001) qtf(j) = summ
end do
!
! compute the qr factorization of the updated jacobian.
!
call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
call r1mpyq(1,n,qtf,1,wa2,wa3)
!
! end of the inner loop.
!
jeval = .false.
end do
!
! end of the outer loop.
!
end do
300 continue
!
! termination, either normal or user imposed.
!
if (iflag < 0) info = iflag
iflag = 0
if (nprint > 0) call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
end subroutine hybrjmv
subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
use const_and_precisions, only : zero, one, epsmch=>comp_eps
! arguments