diff options
Diffstat (limited to 'math/minpack/lmstr.f')
| -rw-r--r-- | math/minpack/lmstr.f | 466 |
1 files changed, 466 insertions, 0 deletions
diff --git a/math/minpack/lmstr.f b/math/minpack/lmstr.f new file mode 100644 index 00000000..d9a7893f --- /dev/null +++ b/math/minpack/lmstr.f @@ -0,0 +1,466 @@ + subroutine lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol, + * maxfev,diag,mode,factor,nprint,info,nfev,njev, + * ipvt,qtf,wa1,wa2,wa3,wa4) + integer m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev + integer ipvt(n) + logical sing + double precision ftol,xtol,gtol,factor + double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),qtf(n), + * wa1(n),wa2(n),wa3(n),wa4(m) +c ********** +c +c subroutine lmstr +c +c the purpose of lmstr is to minimize the sum of the squares of +c m nonlinear functions in n variables by a modification of +c the levenberg-marquardt algorithm which uses minimal storage. +c the user must provide a subroutine which calculates the +c functions and the rows of the jacobian. +c +c the subroutine statement is +c +c subroutine lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol, +c maxfev,diag,mode,factor,nprint,info,nfev, +c njev,ipvt,qtf,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 rows of the jacobian. +c fcn must be declared in an external statement in the +c user calling program, and should be written as follows. +c +c subroutine fcn(m,n,x,fvec,fjrow,iflag) +c integer m,n,iflag +c double precision x(n),fvec(m),fjrow(n) +c ---------- +c if iflag = 1 calculate the functions at x and +c return this vector in fvec. +c if iflag = i calculate the (i-1)-st row of the +c jacobian at x and return this vector in fjrow. +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 lmstr. +c in this case set iflag to a negative integer. +c +c m is a positive integer input variable set to the number +c of functions. +c +c n is a positive integer input variable set to the number +c of variables. n must not exceed m. +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 m which contains +c the functions evaluated at the output x. +c +c fjac is an output n by n array. the upper triangle of fjac +c contains an upper triangular matrix r such that +c +c t t t +c p *(jac *jac)*p = r *r, +c +c where p is a permutation matrix and jac is the final +c calculated jacobian. column j of p is column ipvt(j) +c (see below) of the identity matrix. the lower triangular +c part of fjac contains information generated during +c the computation of r. +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 ftol is a nonnegative input variable. termination +c occurs when both the actual and predicted relative +c reductions in the sum of squares are at most ftol. +c therefore, ftol measures the relative error desired +c in the sum of squares. +c +c xtol is a nonnegative input variable. termination +c occurs when the relative error between two consecutive +c iterates is at most xtol. therefore, xtol measures the +c relative error desired in the approximate solution. +c +c gtol is a nonnegative input variable. termination +c occurs when the cosine of the angle between fvec and +c any column of the jacobian is at most gtol in absolute +c value. therefore, gtol measures the orthogonality +c desired between the function vector and the columns +c of the jacobian. +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. if nprint is not positive, no special calls +c of fcn 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 both actual and predicted relative reductions +c in the sum of squares are at most ftol. +c +c info = 2 relative error between two consecutive iterates +c is at most xtol. +c +c info = 3 conditions for info = 1 and info = 2 both hold. +c +c info = 4 the cosine of the angle between fvec and any +c column of the jacobian is at most gtol in +c absolute value. +c +c info = 5 number of calls to fcn with iflag = 1 has +c reached maxfev. +c +c info = 6 ftol is too small. no further reduction in +c the sum of squares is possible. +c +c info = 7 xtol is too small. no further improvement in +c the approximate solution x is possible. +c +c info = 8 gtol is too small. fvec is orthogonal to the +c columns of the jacobian to machine precision. +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 ipvt is an integer output array of length n. ipvt +c defines a permutation matrix p such that jac*p = q*r, +c where jac is the final calculated jacobian, q is +c orthogonal (not stored), and r is upper triangular. +c column j of p is column ipvt(j) of the identity matrix. +c +c qtf is an output array of length n which contains +c the first n elements of the vector (q transpose)*fvec. +c +c wa1, wa2, and wa3 are work arrays of length n. +c +c wa4 is a work array of length m. +c +c subprograms called +c +c user-supplied ...... fcn +c +c minpack-supplied ... dpmpar,enorm,lmpar,qrfac,rwupdt +c +c fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod +c +c argonne national laboratory. minpack project. march 1980. +c burton s. garbow, dudley v. goetschel, kenneth e. hillstrom, +c jorge j. more +c +c ********** + integer i,iflag,iter,j,l + double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm, + * one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio, + * sum,temp,temp1,temp2,xnorm,zero + double precision dpmpar,enorm + data one,p1,p5,p25,p75,p0001,zero + * /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,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. m .lt. n .or. ldfjac .lt. n + * .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero + * .or. maxfev .le. 0 .or. factor .le. zero) go to 340 + if (mode .ne. 2) go to 20 + do 10 j = 1, n + if (diag(j) .le. zero) go to 340 + 10 continue + 20 continue +c +c evaluate the function at the starting point +c and calculate its norm. +c + iflag = 1 + call fcn(m,n,x,fvec,wa3,iflag) + nfev = 1 + if (iflag .lt. 0) go to 340 + fnorm = enorm(m,fvec) +c +c initialize levenberg-marquardt parameter and iteration counter. +c + par = zero + iter = 1 +c +c beginning of the outer loop. +c + 30 continue +c +c if requested, call fcn to enable printing of iterates. +c + if (nprint .le. 0) go to 40 + iflag = 0 + if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,wa3,iflag) + if (iflag .lt. 0) go to 340 + 40 continue +c +c compute the qr factorization of the jacobian matrix +c calculated one row at a time, while simultaneously +c forming (q transpose)*fvec and storing the first +c n components in qtf. +c + do 60 j = 1, n + qtf(j) = zero + do 50 i = 1, n + fjac(i,j) = zero + 50 continue + 60 continue + iflag = 2 + do 70 i = 1, m + call fcn(m,n,x,fvec,wa3,iflag) + if (iflag .lt. 0) go to 340 + temp = fvec(i) + call rwupdt(n,fjac,ldfjac,wa3,qtf,temp,wa1,wa2) + iflag = iflag + 1 + 70 continue + njev = njev + 1 +c +c if the jacobian is rank deficient, call qrfac to +c reorder its columns and update the components of qtf. +c + sing = .false. + do 80 j = 1, n + if (fjac(j,j) .eq. zero) sing = .true. + ipvt(j) = j + wa2(j) = enorm(j,fjac(1,j)) + 80 continue + if (.not.sing) go to 130 + call qrfac(n,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3) + 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 + fjac(j,j) = wa1(j) + 120 continue + 130 continue +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 170 + if (mode .eq. 2) go to 150 + do 140 j = 1, n + diag(j) = wa2(j) + if (wa2(j) .eq. zero) diag(j) = one + 140 continue + 150 continue +c +c on the first iteration, calculate the norm of the scaled x +c and initialize the step bound delta. +c + do 160 j = 1, n + wa3(j) = diag(j)*x(j) + 160 continue + xnorm = enorm(n,wa3) + delta = factor*xnorm + if (delta .eq. zero) delta = factor + 170 continue +c +c compute the norm of the scaled gradient. +c + gnorm = zero + if (fnorm .eq. zero) go to 210 + do 200 j = 1, n + l = ipvt(j) + if (wa2(l) .eq. zero) go to 190 + sum = zero + do 180 i = 1, j + sum = sum + fjac(i,j)*(qtf(i)/fnorm) + 180 continue + gnorm = dmax1(gnorm,dabs(sum/wa2(l))) + 190 continue + 200 continue + 210 continue +c +c test for convergence of the gradient norm. +c + if (gnorm .le. gtol) info = 4 + if (info .ne. 0) go to 340 +c +c rescale if necessary. +c + if (mode .eq. 2) go to 230 + do 220 j = 1, n + diag(j) = dmax1(diag(j),wa2(j)) + 220 continue + 230 continue +c +c beginning of the inner loop. +c + 240 continue +c +c determine the levenberg-marquardt parameter. +c + call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2, + * wa3,wa4) +c +c store the direction p and x + p. calculate the norm of p. +c + do 250 j = 1, n + wa1(j) = -wa1(j) + wa2(j) = x(j) + wa1(j) + wa3(j) = diag(j)*wa1(j) + 250 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(m,n,wa2,wa4,wa3,iflag) + nfev = nfev + 1 + if (iflag .lt. 0) go to 340 + fnorm1 = enorm(m,wa4) +c +c compute the scaled actual reduction. +c + actred = -one + if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2 +c +c compute the scaled predicted reduction and +c the scaled directional derivative. +c + do 270 j = 1, n + wa3(j) = zero + l = ipvt(j) + temp = wa1(l) + do 260 i = 1, j + wa3(i) = wa3(i) + fjac(i,j)*temp + 260 continue + 270 continue + temp1 = enorm(n,wa3)/fnorm + temp2 = (dsqrt(par)*pnorm)/fnorm + prered = temp1**2 + temp2**2/p5 + dirder = -(temp1**2 + temp2**2) +c +c compute the ratio of the actual to the predicted +c reduction. +c + ratio = zero + if (prered .ne. zero) ratio = actred/prered +c +c update the step bound. +c + if (ratio .gt. p25) go to 280 + if (actred .ge. zero) temp = p5 + if (actred .lt. zero) + * temp = p5*dirder/(dirder + p5*actred) + if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1 + delta = temp*dmin1(delta,pnorm/p1) + par = par/temp + go to 300 + 280 continue + if (par .ne. zero .and. ratio .lt. p75) go to 290 + delta = pnorm/p5 + par = p5*par + 290 continue + 300 continue +c +c test for successful iteration. +c + if (ratio .lt. p0001) go to 330 +c +c successful iteration. update x, fvec, and their norms. +c + do 310 j = 1, n + x(j) = wa2(j) + wa2(j) = diag(j)*x(j) + 310 continue + do 320 i = 1, m + fvec(i) = wa4(i) + 320 continue + xnorm = enorm(n,wa2) + fnorm = fnorm1 + iter = iter + 1 + 330 continue +c +c tests for convergence. +c + if (dabs(actred) .le. ftol .and. prered .le. ftol + * .and. p5*ratio .le. one) info = 1 + if (delta .le. xtol*xnorm) info = 2 + if (dabs(actred) .le. ftol .and. prered .le. ftol + * .and. p5*ratio .le. one .and. info .eq. 2) info = 3 + if (info .ne. 0) go to 340 +c +c tests for termination and stringent tolerances. +c + if (nfev .ge. maxfev) info = 5 + if (dabs(actred) .le. epsmch .and. prered .le. epsmch + * .and. p5*ratio .le. one) info = 6 + if (delta .le. epsmch*xnorm) info = 7 + if (gnorm .le. epsmch) info = 8 + if (info .ne. 0) go to 340 +c +c end of the inner loop. repeat if iteration unsuccessful. +c + if (ratio .lt. p0001) go to 240 +c +c end of the outer loop. +c + go to 30 + 340 continue +c +c termination, either normal or user imposed. +c + if (iflag .lt. 0) info = iflag + iflag = 0 + if (nprint .gt. 0) call fcn(m,n,x,fvec,wa3,iflag) + return +c +c last card of subroutine lmstr. +c + end |