hybrd – fortran-astrodynamics-toolkit

public subroutine hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4)

The purpose of hybrd 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. the jacobian is then calculated by a forward-difference approximation.

Characteristics of the algorithm.

HYBRD is a modification of the Powell hybrid method. Two of its main characteristics involve the choice of the correction as a convex combination of the Newton and scaled gradient directions and the updating of the Jacobian by the rank-1 method of Broy- den. The choice of the correction guarantees (under reasonable conditions) global convergence for starting points far from the solution and a fast rate of convergence. The Jacobian is approximated by forward differences at the starting point, but forward differences are not used again until the rank-1 method fails to produce satisfactory progress.

References

M. J. D. Powell, A Hybrid Method for Nonlinear Equations. Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970.

Arguments

Type	Intent		Name
procedure(fcn_hybrd)		::	fcn	user-supplied subroutine which calculates the functions
integer,	intent(in)	::	n	a positive integer input variable set to the number of functions and variables.
real(kind=wp),	intent(inout)	::	x(n)	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.
real(kind=wp),	intent(out)	::	fvec(n)	an output array of length `n` which contains the functions evaluated at the output `x`.
real(kind=wp),	intent(in)	::	xtol	a nonnegative input variable. termination occurs when the relative error between two consecutive iterates is at most `xtol`.
integer,	intent(in)	::	maxfev	a positive integer input variable. termination occurs when the number of calls to `fcn` is at least `maxfev` by the end of an iteration.
integer,	intent(in)	::	ml	a nonnegative integer input variable which specifies the number of subdiagonals within the band of the jacobian matrix. if the jacobian is not banded, set `ml` to at least `n - 1`.
integer,	intent(in)	::	mu	a nonnegative integer input variable which specifies the number of superdiagonals within the band of the jacobian matrix. if the jacobian is not banded, set `mu` to at least`n - 1`.
real(kind=wp),	intent(in)	::	epsfcn	an input variable used in determining a suitable step length for the forward-difference approximation. this approximation assumes that the relative errors in the functions are of the order of `epsfcn`. if `epsfcn` is less than the machine precision, it is assumed that the relative errors in the functions are of the order of the machine precision.
real(kind=wp),	intent(inout)	::	diag(n)	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.
integer,	intent(in)	::	mode	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`.
real(kind=wp),	intent(in)	::	factor	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.
integer,	intent(in)	::	nprint	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. if `nprint` is not positive, no special calls of `fcn` with `iflag = 0` are made.
integer,	intent(out)	::	info	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` has reached or exceeded `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.
integer,	intent(out)	::	nfev	output variable set to the number of calls to `fcn`.
real(kind=wp),	intent(out)	::	fjac(ldfjac,n)	array which contains the orthogonal matrix `q` produced by the QR factorization of the final approximate jacobian.
integer,	intent(in)	::	ldfjac	a positive integer input variable not less than `n` which specifies the leading dimension of the array `fjac`.
real(kind=wp),	intent(out)	::	r(lr)	an output array which contains the upper triangular matrix produced by the QR factorization of the final approximate jacobian, stored rowwise.
integer,	intent(in)	::	lr	a positive integer input variable not less than `(n*(n+1))/2`.
real(kind=wp),	intent(out)	::	qtf(n)	an output array of length `n` which contains the vector `(q transpose)*fvec`.
real(kind=wp),	intent(inout)	::	wa1(n)	work array
real(kind=wp),	intent(inout)	::	wa2(n)	work array
real(kind=wp),	intent(inout)	::	wa3(n)	work array
real(kind=wp),	intent(inout)	::	wa4(n)	work array

Calls

Help

Called by

Help

Source Code

    subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,mode, &
                     factor,nprint,info,nfev,fjac,ldfjac,r,lr,qtf,wa1,&
                     wa2,wa3,wa4)

    implicit none

    procedure(fcn_hybrd) :: fcn             !! user-supplied subroutine which calculates the functions
    integer,intent(in) :: n                 !! a positive integer input variable set to the number
                                            !! of functions and variables.
    integer,intent(in) :: maxfev            !! a positive integer input variable. termination
                                            !! occurs when the number of calls to `fcn` is at least `maxfev`
                                            !! by the end of an iteration.
    integer,intent(in) :: ml                !! a nonnegative integer input variable which specifies
                                            !! the number of subdiagonals within the band of the
                                            !! jacobian matrix. if the jacobian is not banded, set
                                            !! `ml` to at least `n - 1`.
    integer,intent(in) :: mu                !! a nonnegative integer input variable which specifies
                                            !! the number of superdiagonals within the band of the
                                            !! jacobian matrix. if the jacobian is not banded, set
                                            !! `mu` to at least` n - 1`.
    integer,intent(in) :: mode              !! 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`.
    integer,intent(in)  :: nprint           !! 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. if `nprint` is not positive, no special calls
                                            !! of `fcn` with `iflag = 0` are made.
    integer,intent(out) :: info             !! 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` has reached or exceeded
                                            !!    `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.
    integer,intent(out) :: nfev             !! output variable set to the number of calls to `fcn`.
    integer,intent(in):: ldfjac             !! a positive integer input variable not less than `n`
                                            !! which specifies the leading dimension of the array `fjac`.
    integer,intent(in) :: lr                !! a positive integer input variable not less than `(n*(n+1))/2`.
    real(wp),intent(in) :: xtol             !! a nonnegative input variable. termination
                                            !! occurs when the relative error between two consecutive
                                            !! iterates is at most `xtol`.
    real(wp),intent(in) :: epsfcn           !! an input variable used in determining a suitable
                                            !! step length for the forward-difference approximation. this
                                            !! approximation assumes that the relative errors in the
                                            !! functions are of the order of `epsfcn`. if `epsfcn` is less
                                            !! than the machine precision, it is assumed that the relative
                                            !! errors in the functions are of the order of the machine
                                            !! precision.
    real(wp),intent(in) :: factor           !! 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.
    real(wp),intent(inout) :: x(n)          !! 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.
    real(wp),intent(out) :: fvec(n)         !! an output array of length `n` which contains
                                            !! the functions evaluated at the output `x`.
    real(wp),intent(inout) :: diag(n)       !! 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.
    real(wp),intent(out) :: fjac(ldfjac,n)  !! array which contains the
                                            !! orthogonal matrix `q` produced by the QR factorization
                                            !! of the final approximate jacobian.
    real(wp),intent(out) :: r(lr)           !! an output array which contains the
                                            !! upper triangular matrix produced by the QR factorization
                                            !! of the final approximate jacobian, stored rowwise.
    real(wp),intent(out) :: qtf(n)          !! an output array of length `n` which contains
                                            !! the vector `(q transpose)*fvec`.
    real(wp),intent(inout) :: wa1(n)  !! work array
    real(wp),intent(inout) :: wa2(n)  !! work array
    real(wp),intent(inout) :: wa3(n)  !! work array
    real(wp),intent(inout) :: wa4(n)  !! work array

    integer :: i , iflag , iter , j , jm1 , l , msum , ncfail , ncsuc , nslow1 , nslow2
    integer :: iwa(1)
    logical :: jeval , sing
    real(wp) :: actred , delta , epsmch , fnorm , fnorm1 , &
                  pnorm , prered , ratio ,&
                  sum , temp , xnorm

    real(wp),parameter :: p1    = 1.0e-1_wp
    real(wp),parameter :: p5    = 5.0e-1_wp
    real(wp),parameter :: p001  = 1.0e-3_wp
    real(wp),parameter :: p0001 = 1.0e-4_wp

    epsmch = dpmpar(1)  ! the machine precision

    info = 0
    iflag = 0
    nfev = 0
    !
    !     check the input parameters for errors.
    !
    if ( n<=0 .or. xtol<zero .or. maxfev<=0 .or. ml<0 .or. mu<0 .or.  &
         factor<=zero .or. ldfjac<n .or. lr<(n*(n+1))/2 ) goto 300
    if ( mode==2 ) then
       do j = 1 , n
          if ( diag(j)<=zero ) goto 300
       enddo
    endif
    !
    !     evaluate the function at the starting point
    !     and calculate its norm.
    !
    iflag = 1
    call fcn(n,x,fvec,iflag)
    nfev = 1
    if ( iflag<0 ) goto 300
    fnorm = enorm(n,fvec)
    !
    !     determine the number of calls to fcn needed to compute
    !     the jacobian matrix.
    !
    msum = min0(ml+mu+1,n)
    !
    !     initialize iteration counter and monitors.
    !
    iter = 1
    ncsuc = 0
    ncfail = 0
    nslow1 = 0
    nslow2 = 0
    !
    !     beginning of the outer loop.
    !
    100  jeval = .true.
    !
    !        calculate the jacobian matrix.
    !
    iflag = 2
    call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,wa2)
    nfev = nfev + msum
    if ( iflag<0 ) goto 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
          enddo
       endif
    !
    !        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)
       enddo
       xnorm = enorm(n,wa3)
       delta = factor*xnorm
       if ( delta==zero ) delta = factor
    endif
    !
    !        form (q transpose)*fvec and store in qtf.
    !
    do i = 1 , n
       qtf(i) = fvec(i)
    enddo
    do j = 1 , n
       if ( fjac(j,j)/=zero ) then
          sum = zero
          do i = j , n
             sum = sum + fjac(i,j)*qtf(i)
          enddo
          temp = -sum/fjac(j,j)
          do i = j , n
             qtf(i) = qtf(i) + fjac(i,j)*temp
          enddo
       endif
    enddo
    !
    !        copy the triangular factor of the qr factorization into r.
    !
    sing = .false.
    do j = 1 , n
       l = j
       jm1 = j - 1
       if ( jm1>=1 ) then
          do i = 1 , jm1
             r(l) = fjac(i,j)
             l = l + n - i
          enddo
       endif
       r(l) = wa1(j)
       if ( wa1(j)==zero ) sing = .true.
    enddo
    !
    !        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))
       enddo
    endif
    !
    !        beginning of the inner loop.
    !
    !
    !           if requested, call fcn to enable printing of iterates.
    !
    200  if ( nprint>0 ) then
       iflag = 0
       if ( mod(iter-1,nprint)==0 ) call fcn(n,x,fvec,iflag)
       if ( iflag<0 ) goto 300
    endif
    !
    !           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)
    enddo
    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,wa4,iflag)
    nfev = nfev + 1
    if ( iflag>=0 ) then
       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
          sum = zero
          do j = i , n
             sum = sum + r(l)*wa1(j)
             l = l + 1
          enddo
          wa3(i) = qtf(i) + sum
       enddo
       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
          ncfail = 0
          ncsuc = ncsuc + 1
          if ( ratio>=p5 .or. ncsuc>1 ) delta = dmax1(delta,pnorm/p5)
          if ( abs(ratio-one)<=p1 ) delta = pnorm/p5
       else
          ncsuc = 0
          ncfail = ncfail + 1
          delta = p5*delta
       endif
    !
    !           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)
          enddo
          xnorm = enorm(n,wa2)
          fnorm = fnorm1
          iter = iter + 1
       endif
    !
    !           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 ) then
    !
    !           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 ) then
    !
    !           criterion for recalculating jacobian approximation
    !           by forward differences.
    !
             if ( ncfail==2 ) goto 100
    !
    !           calculate the rank one modification to the jacobian
    !           and update qtf if necessary.
    !
             do j = 1 , n
                sum = zero
                do i = 1 , n
                   sum = sum + fjac(i,j)*wa4(i)
                enddo
                wa2(j) = (sum-wa3(j))/pnorm
                wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
                if ( ratio>=p0001 ) qtf(j) = sum
             enddo
    !
    !           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 of the outer loop.
    !
             goto 200
          endif
       endif
    endif
    !
    !     termination, either normal or user imposed.
    !
    300  if ( iflag<0 ) info = iflag
    iflag = 0
    if ( nprint>0 ) call fcn(n,x,fvec,iflag)

    end subroutine hybrd

hybrd Subroutine

Contents

Source Code

public subroutine hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4)

Characteristics of the algorithm.

References

Arguments

Calls

Called by

Source Code