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 real xtol,factor real 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 real 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,spmpar,enorm, c qform,qrfac,r1mpyq,r1updt c c fortran-supplied ... abs,amax1,amin1,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 real actred,delta,epsmch,fnorm,fnorm1,one,pnorm,prered,p1,p5, * p001,p0001,ratio,sum,temp,xnorm,zero real spmpar,enorm data one,p1,p5,p001,p0001,zero * /1.0e0,1.0e-1,5.0e-1,1.0e-3,1.0e-4,0.0e0/ c c epsmch is the machine precision. c epsmch = spmpar(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) = amax1(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 = amin1(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 = amax1(delta,pnorm/p5) if (abs(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*amax1(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