diff options
Diffstat (limited to 'math/minpack/hybrj.f')
| -rw-r--r-- | math/minpack/hybrj.f | 440 |
1 files changed, 440 insertions, 0 deletions
diff --git a/math/minpack/hybrj.f b/math/minpack/hybrj.f new file mode 100644 index 00000000..3070dad3 --- /dev/null +++ b/math/minpack/hybrj.f @@ -0,0 +1,440 @@ + subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, + * factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, + * wa3,wa4) + integer n,ldfjac,maxfev,mode,nprint,info,nfev,njev,lr + double precision xtol,factor + double precision x(n),fvec(n),fjac(ldfjac,n),diag(n),r(lr), + * qtf(n),wa1(n),wa2(n),wa3(n),wa4(n) +c ********** +c +c subroutine hybrj +c +c the purpose of hybrj is to find a zero of a system of +c n nonlinear functions in n variables by a modification +c of the powell hybrid method. the user must provide a +c subroutine which calculates the functions and the jacobian. +c +c the subroutine statement is +c +c subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag, +c mode,factor,nprint,info,nfev,njev,r,lr,qtf, +c wa1,wa2,wa3,wa4) +c +c where +c +c fcn is the name of the user-supplied subroutine which +c calculates the functions and the jacobian. fcn must +c be declared in an external statement in the user +c calling program, and should be written as follows. +c +c subroutine fcn(n,x,fvec,fjac,ldfjac,iflag) +c integer n,ldfjac,iflag +c double precision x(n),fvec(n),fjac(ldfjac,n) +c ---------- +c if iflag = 1 calculate the functions at x and +c return this vector in fvec. do not alter fjac. +c if iflag = 2 calculate the jacobian at x and +c return this matrix in fjac. do not alter fvec. +c --------- +c return +c end +c +c the value of iflag should not be changed by fcn unless +c the user wants to terminate execution of hybrj. +c in this case set iflag to a negative integer. +c +c n is a positive integer input variable set to the number +c of functions and variables. +c +c x is an array of length n. on input x must contain +c an initial estimate of the solution vector. on output x +c contains the final estimate of the solution vector. +c +c fvec is an output array of length n which contains +c the functions evaluated at the output x. +c +c fjac is an output n by n array which contains the +c orthogonal matrix q produced by the qr factorization +c of the final approximate jacobian. +c +c ldfjac is a positive integer input variable not less than n +c which specifies the leading dimension of the array fjac. +c +c xtol is a nonnegative input variable. termination +c occurs when the relative error between two consecutive +c iterates is at most xtol. +c +c maxfev is a positive integer input variable. termination +c occurs when the number of calls to fcn with iflag = 1 +c has reached maxfev. +c +c diag is an array of length n. if mode = 1 (see +c below), diag is internally set. if mode = 2, diag +c must contain positive entries that serve as +c multiplicative scale factors for the variables. +c +c mode is an integer input variable. if mode = 1, the +c variables will be scaled internally. if mode = 2, +c the scaling is specified by the input diag. other +c values of mode are equivalent to mode = 1. +c +c factor is a positive input variable used in determining the +c initial step bound. this bound is set to the product of +c factor and the euclidean norm of diag*x if nonzero, or else +c to factor itself. in most cases factor should lie in the +c interval (.1,100.). 100. is a generally recommended value. +c +c nprint is an integer input variable that enables controlled +c printing of iterates if it is positive. in this case, +c fcn is called with iflag = 0 at the beginning of the first +c iteration and every nprint iterations thereafter and +c immediately prior to return, with x and fvec available +c for printing. fvec and fjac should not be altered. +c if nprint is not positive, no special calls of fcn +c with iflag = 0 are made. +c +c info is an integer output variable. if the user has +c terminated execution, info is set to the (negative) +c value of iflag. see description of fcn. otherwise, +c info is set as follows. +c +c info = 0 improper input parameters. +c +c info = 1 relative error between two consecutive iterates +c is at most xtol. +c +c info = 2 number of calls to fcn with iflag = 1 has +c reached maxfev. +c +c info = 3 xtol is too small. no further improvement in +c the approximate solution x is possible. +c +c info = 4 iteration is not making good progress, as +c measured by the improvement from the last +c five jacobian evaluations. +c +c info = 5 iteration is not making good progress, as +c measured by the improvement from the last +c ten iterations. +c +c nfev is an integer output variable set to the number of +c calls to fcn with iflag = 1. +c +c njev is an integer output variable set to the number of +c calls to fcn with iflag = 2. +c +c r is an output array of length lr which contains the +c upper triangular matrix produced by the qr factorization +c of the final approximate jacobian, stored rowwise. +c +c lr is a positive integer input variable not less than +c (n*(n+1))/2. +c +c qtf is an output array of length n which contains +c the vector (q transpose)*fvec. +c +c wa1, wa2, wa3, and wa4 are work arrays of length n. +c +c subprograms called +c +c user-supplied ...... fcn +c +c minpack-supplied ... dogleg,dpmpar,enorm, +c qform,qrfac,r1mpyq,r1updt +c +c fortran-supplied ... dabs,dmax1,dmin1,mod +c +c argonne national laboratory. minpack project. march 1980. +c burton s. garbow, kenneth e. hillstrom, jorge j. more +c +c ********** + integer i,iflag,iter,j,jm1,l,ncfail,ncsuc,nslow1,nslow2 + integer iwa(1) + logical jeval,sing + double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm, + * prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm, + * zero + double precision dpmpar,enorm + data one,p1,p5,p001,p0001,zero + * /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/ +c +c epsmch is the machine precision. +c + epsmch = dpmpar(1) +c + info = 0 + iflag = 0 + nfev = 0 + njev = 0 +c +c check the input parameters for errors. +c + if (n .le. 0 .or. ldfjac .lt. n .or. xtol .lt. zero + * .or. maxfev .le. 0 .or. factor .le. zero + * .or. lr .lt. (n*(n + 1))/2) go to 300 + if (mode .ne. 2) go to 20 + do 10 j = 1, n + if (diag(j) .le. zero) go to 300 + 10 continue + 20 continue +c +c evaluate the function at the starting point +c and calculate its norm. +c + iflag = 1 + call fcn(n,x,fvec,fjac,ldfjac,iflag) + nfev = 1 + if (iflag .lt. 0) go to 300 + fnorm = enorm(n,fvec) +c +c initialize iteration counter and monitors. +c + iter = 1 + ncsuc = 0 + ncfail = 0 + nslow1 = 0 + nslow2 = 0 +c +c beginning of the outer loop. +c + 30 continue + jeval = .true. +c +c calculate the jacobian matrix. +c + iflag = 2 + call fcn(n,x,fvec,fjac,ldfjac,iflag) + njev = njev + 1 + if (iflag .lt. 0) go to 300 +c +c compute the qr factorization of the jacobian. +c + call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3) +c +c on the first iteration and if mode is 1, scale according +c to the norms of the columns of the initial jacobian. +c + if (iter .ne. 1) go to 70 + if (mode .eq. 2) go to 50 + do 40 j = 1, n + diag(j) = wa2(j) + if (wa2(j) .eq. zero) diag(j) = one + 40 continue + 50 continue +c +c on the first iteration, calculate the norm of the scaled x +c and initialize the step bound delta. +c + do 60 j = 1, n + wa3(j) = diag(j)*x(j) + 60 continue + xnorm = enorm(n,wa3) + delta = factor*xnorm + if (delta .eq. zero) delta = factor + 70 continue +c +c form (q transpose)*fvec and store in qtf. +c + do 80 i = 1, n + qtf(i) = fvec(i) + 80 continue + do 120 j = 1, n + if (fjac(j,j) .eq. zero) go to 110 + sum = zero + do 90 i = j, n + sum = sum + fjac(i,j)*qtf(i) + 90 continue + temp = -sum/fjac(j,j) + do 100 i = j, n + qtf(i) = qtf(i) + fjac(i,j)*temp + 100 continue + 110 continue + 120 continue +c +c copy the triangular factor of the qr factorization into r. +c + sing = .false. + do 150 j = 1, n + l = j + jm1 = j - 1 + if (jm1 .lt. 1) go to 140 + do 130 i = 1, jm1 + r(l) = fjac(i,j) + l = l + n - i + 130 continue + 140 continue + r(l) = wa1(j) + if (wa1(j) .eq. zero) sing = .true. + 150 continue +c +c accumulate the orthogonal factor in fjac. +c + call qform(n,n,fjac,ldfjac,wa1) +c +c rescale if necessary. +c + if (mode .eq. 2) go to 170 + do 160 j = 1, n + diag(j) = dmax1(diag(j),wa2(j)) + 160 continue + 170 continue +c +c beginning of the inner loop. +c + 180 continue +c +c if requested, call fcn to enable printing of iterates. +c + if (nprint .le. 0) go to 190 + iflag = 0 + if (mod(iter-1,nprint) .eq. 0) + * call fcn(n,x,fvec,fjac,ldfjac,iflag) + if (iflag .lt. 0) go to 300 + 190 continue +c +c determine the direction p. +c + call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3) +c +c store the direction p and x + p. calculate the norm of p. +c + do 200 j = 1, n + wa1(j) = -wa1(j) + wa2(j) = x(j) + wa1(j) + wa3(j) = diag(j)*wa1(j) + 200 continue + pnorm = enorm(n,wa3) +c +c on the first iteration, adjust the initial step bound. +c + if (iter .eq. 1) delta = dmin1(delta,pnorm) +c +c evaluate the function at x + p and calculate its norm. +c + iflag = 1 + call fcn(n,wa2,wa4,fjac,ldfjac,iflag) + nfev = nfev + 1 + if (iflag .lt. 0) go to 300 + fnorm1 = enorm(n,wa4) +c +c compute the scaled actual reduction. +c + actred = -one + if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2 +c +c compute the scaled predicted reduction. +c + l = 1 + do 220 i = 1, n + sum = zero + do 210 j = i, n + sum = sum + r(l)*wa1(j) + l = l + 1 + 210 continue + wa3(i) = qtf(i) + sum + 220 continue + temp = enorm(n,wa3) + prered = zero + if (temp .lt. fnorm) prered = one - (temp/fnorm)**2 +c +c compute the ratio of the actual to the predicted +c reduction. +c + ratio = zero + if (prered .gt. zero) ratio = actred/prered +c +c update the step bound. +c + if (ratio .ge. p1) go to 230 + ncsuc = 0 + ncfail = ncfail + 1 + delta = p5*delta + go to 240 + 230 continue + ncfail = 0 + ncsuc = ncsuc + 1 + if (ratio .ge. p5 .or. ncsuc .gt. 1) + * delta = dmax1(delta,pnorm/p5) + if (dabs(ratio-one) .le. p1) delta = pnorm/p5 + 240 continue +c +c test for successful iteration. +c + if (ratio .lt. p0001) go to 260 +c +c successful iteration. update x, fvec, and their norms. +c + do 250 j = 1, n + x(j) = wa2(j) + wa2(j) = diag(j)*x(j) + fvec(j) = wa4(j) + 250 continue + xnorm = enorm(n,wa2) + fnorm = fnorm1 + iter = iter + 1 + 260 continue +c +c determine the progress of the iteration. +c + nslow1 = nslow1 + 1 + if (actred .ge. p001) nslow1 = 0 + if (jeval) nslow2 = nslow2 + 1 + if (actred .ge. p1) nslow2 = 0 +c +c test for convergence. +c + if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1 + if (info .ne. 0) go to 300 +c +c tests for termination and stringent tolerances. +c + if (nfev .ge. maxfev) info = 2 + if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3 + if (nslow2 .eq. 5) info = 4 + if (nslow1 .eq. 10) info = 5 + if (info .ne. 0) go to 300 +c +c criterion for recalculating jacobian. +c + if (ncfail .eq. 2) go to 290 +c +c calculate the rank one modification to the jacobian +c and update qtf if necessary. +c + do 280 j = 1, n + sum = zero + do 270 i = 1, n + sum = sum + fjac(i,j)*wa4(i) + 270 continue + wa2(j) = (sum - wa3(j))/pnorm + wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm) + if (ratio .ge. p0001) qtf(j) = sum + 280 continue +c +c compute the qr factorization of the updated jacobian. +c + 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) +c +c end of the inner loop. +c + jeval = .false. + go to 180 + 290 continue +c +c end of the outer loop. +c + go to 30 + 300 continue +c +c termination, either normal or user imposed. +c + if (iflag .lt. 0) info = iflag + iflag = 0 + if (nprint .gt. 0) call fcn(n,x,fvec,fjac,ldfjac,iflag) + return +c +c last card of subroutine hybrj. +c + end |