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.
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,fvec,fjac,ldfjac,iflag)
integer n,ldfjac,iflag
real(wp) 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 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.
HYBRJ 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 calcu- lated at the starting point, but it is not recalculated until the rank-1 method fails to produce satisfactory progress.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| procedure(fcn_hybrj) | :: | fcn | ||||
| integer | :: | n | ||||
| real(kind=wp) | :: | x(n) | ||||
| real(kind=wp) | :: | fvec(n) | ||||
| real(kind=wp) | :: | fjac(ldfjac,n) | ||||
| integer | :: | ldfjac | ||||
| real(kind=wp) | :: | xtol | ||||
| integer | :: | maxfev | ||||
| real(kind=wp) | :: | diag(n) | ||||
| integer | :: | mode | ||||
| real(kind=wp) | :: | factor | ||||
| integer | :: | nprint | ||||
| integer | :: | info | ||||
| integer | :: | nfev | ||||
| integer | :: | njev | ||||
| real(kind=wp) | :: | r(lr) | ||||
| integer | :: | lr | ||||
| real(kind=wp) | :: | qtf(n) | ||||
| real(kind=wp) | :: | wa1(n) | ||||
| real(kind=wp) | :: | wa2(n) | ||||
| real(kind=wp) | :: | wa3(n) | ||||
| real(kind=wp) | :: | wa4(n) |
subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, & & factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, & & wa3,wa4) implicit none procedure(fcn_hybrj) :: fcn integer n , ldfjac , maxfev , mode , nprint , info , nfev , njev , lr real(wp) xtol , factor real(wp) x(n) , fvec(n) , fjac(ldfjac,n) , diag(n) , & r(lr) , qtf(n) , wa1(n) , wa2(n) , wa3(n) , & wa4(n) integer i , iflag , iter , j , jm1 , l , 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 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 ) 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,fjac,ldfjac,iflag) nfev = 1 if ( iflag<0 ) goto 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. ! 100 jeval = .true. ! ! calculate the jacobian matrix. ! iflag = 2 call fcn(n,x,fvec,fjac,ldfjac,iflag) njev = njev + 1 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) = max(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,fjac,ldfjac,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 = min(delta,pnorm) ! ! evaluate the function at x + p and calculate its norm. ! iflag = 1 call fcn(n,wa2,wa4,fjac,ldfjac,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 = max(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*max(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. ! 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,fjac,ldfjac,iflag) end subroutine hybrj