From 15a1f866b4384efc0bb4b60ee3111e0c80f1c63c Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 27 Aug 2024 18:19:28 +0200 Subject: [PATCH] src/equilibrium: rewrite points_tgo, points_ox 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. --- src/equilibrium.f90 | 225 ++++++++-------- src/magsurf_data.f90 | 8 +- src/vendor/minpack.f90 | 578 ----------------------------------------- 3 files changed, 108 insertions(+), 703 deletions(-) diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 4e386b9..f26bee7 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -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 diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index 027e54e..58e8a87 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -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 diff --git a/src/vendor/minpack.f90 b/src/vendor/minpack.f90 index a30a800..c230b48 100644 --- a/src/vendor/minpack.f90 +++ b/src/vendor/minpack.f90 @@ -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