!< author: John Burkardt !< !< <span style="color: #337ab7; font-family: cabin; font-size: 1.5em;"> !< **A fortran collection of functions for minimization problems** !< </span> !< !< @note !< MINPACK is a FORTRAN90 library which solves systems of nonlinear equations, !< or carries out the least squares minimization of the residual of a set of linear or nonlinear equations, !< by Jorge More, Danny Sorenson, Burton Garbow, Kenneth Hillstrom. !< !< [link](https://people.math.sc.edu/Burkardt/f_src/minpack/minpack.html) !< !< @endnote !< module minpack implicit none public contains subroutine chkder ( m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err ) !*****************************************************************************80 ! !! CHKDER checks the gradients of M functions of N variables. ! ! Discussion: ! ! CHKDER checks the gradients of M nonlinear functions in N variables, ! evaluated at a point X, for consistency with the functions themselves. ! ! The user calls CHKDER twice, first with MODE = 1 and then with MODE = 2. ! ! MODE = 1. ! On input, ! X contains the point of evaluation. ! On output, ! XP is set to a neighboring point. ! ! Now the user must evaluate the function and gradients at X, and the ! function at XP. Then the function is called again: ! ! MODE = 2. ! On input, ! FVEC contains the function values at X, ! FJAC contains the function gradients at X. ! FVECP contains the functions evaluated at XP. ! On output, ! ERR contains measures of correctness of the respective gradients. ! ! The function does not perform reliably if cancellation or ! rounding errors cause a severe loss of significance in the ! evaluation of a function. Therefore, none of the components ! of X should be unusually small (in particular, zero) or any ! other value which may cause loss of significance. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, is the number of functions. ! ! Input, integer ( kind = 4 ) N, is the number of variables. ! ! Input, real ( kind = 8 ) X(N), the point at which the jacobian is to be ! evaluated. ! ! Input, real ( kind = 8 ) FVEC(M), is used only when MODE = 2. ! In that case, it should contain the function values at X. ! ! Input, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. When MODE = 2, ! FJAC(I,J) should contain the value of dF(I)/dX(J). ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least M. ! ! Output, real ( kind = 8 ) XP(N), on output with MODE = 1, is a neighboring ! point of X, at which the function is to be evaluated. ! ! Input, real ( kind = 8 ) FVECP(M), on input with MODE = 2, is the function ! value at XP. ! ! Input, integer ( kind = 4 ) MODE, should be set to 1 on the first call, and ! 2 on the second. ! ! Output, real ( kind = 8 ) ERR(M). On output when MODE = 2, ERR contains ! measures of correctness of the respective gradients. If there is no ! severe loss of significance, then if ERR(I): ! = 1.0D+00, the I-th gradient is correct, ! = 0.0D+00, the I-th gradient is incorrect. ! > 0.5D+00, the I-th gradient is probably correct. ! < 0.5D+00, the I-th gradient is probably incorrect. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) eps real ( kind = 8 ) epsf real ( kind = 8 ) epslog real ( kind = 8 ) epsmch real ( kind = 8 ) err(m) real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fvec(m) real ( kind = 8 ) fvecp(m) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) mode real ( kind = 8 ) temp real ( kind = 8 ) x(n) real ( kind = 8 ) xp(n) epsmch = epsilon ( epsmch ) eps = sqrt ( epsmch ) ! ! MODE = 1. ! if ( mode == 1 ) then do j = 1, n temp = eps * abs ( x(j) ) if ( temp == 0.0D+00 ) then temp = eps end if xp(j) = x(j) + temp end do ! ! MODE = 2. ! else if ( mode == 2 ) then epsf = 100.0D+00 * epsmch epslog = log10 ( eps ) err = 0.0D+00 do j = 1, n temp = abs ( x(j) ) if ( temp == 0.0D+00 ) then temp = 1.0D+00 end if err(1:m) = err(1:m) + temp * fjac(1:m,j) end do do i = 1, m temp = 1.0D+00 if ( fvec(i) /= 0.0D+00 .and. fvecp(i) /= 0.0D+00 .and. & abs ( fvecp(i)-fvec(i)) >= epsf * abs ( fvec(i) ) ) then temp = eps * abs ( (fvecp(i)-fvec(i)) / eps - err(i) ) & / ( abs ( fvec(i) ) + abs ( fvecp(i) ) ) end if err(i) = 1.0D+00 if ( epsmch < temp .and. temp < eps ) then err(i) = ( log10 ( temp ) - epslog ) / epslog end if if ( eps <= temp ) then err(i) = 0.0D+00 end if end do end if return endsubroutine chkder subroutine dogleg ( n, r, lr, diag, qtb, delta, x ) !*****************************************************************************80 ! !! DOGLEG finds the minimizing combination of Gauss-Newton and gradient steps. ! ! Discussion: ! ! Given an M by N matrix A, an N by N nonsingular diagonal ! matrix D, an M-vector B, and a positive number DELTA, the ! problem is to determine the convex combination X of the ! Gauss-Newton and scaled gradient directions that minimizes ! (A*X - B) in the least squares sense, subject to the ! restriction that the euclidean norm of D*X be at most DELTA. ! ! This function completes the solution of the problem ! if it is provided with the necessary information from the ! QR factorization of A. That is, if A = Q*R, where Q has ! orthogonal columns and R is an upper triangular matrix, ! then DOGLEG expects the full upper triangle of R and ! the first N components of Q'*B. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix R. ! ! Input, real ( kind = 8 ) R(LR), the upper triangular matrix R stored ! by rows. ! ! Input, integer ( kind = 4 ) LR, the size of the R array, which must be ! no less than (N*(N+1))/2. ! ! Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D. ! ! Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'* B. ! ! Input, real ( kind = 8 ) DELTA, is a positive upper bound on the ! euclidean norm of D*X(1:N). ! ! Output, real ( kind = 8 ) X(N), the desired convex combination of the ! Gauss-Newton direction and the scaled gradient direction. ! implicit none integer ( kind = 4 ) lr integer ( kind = 4 ) n real ( kind = 8 ) alpha real ( kind = 8 ) bnorm real ( kind = 8 ) delta real ( kind = 8 ) diag(n) !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsmch real ( kind = 8 ) gnorm integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) k integer ( kind = 4 ) l real ( kind = 8 ) qnorm real ( kind = 8 ) qtb(n) real ( kind = 8 ) r(lr) real ( kind = 8 ) sgnorm real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) x(n) epsmch = epsilon ( epsmch ) ! ! Calculate the Gauss-Newton direction. ! jj = ( n * ( n + 1 ) ) / 2 + 1 do k = 1, n j = n - k + 1 jj = jj - k l = jj + 1 sum2 = 0.0D+00 do i = j + 1, n sum2 = sum2 + r(l) * x(i) l = l + 1 end do temp = r(jj) if ( temp == 0.0D+00 ) then l = j do i = 1, j temp = max ( temp, abs ( r(l)) ) l = l + n - i end do if ( temp == 0.0D+00 ) then temp = epsmch else temp = epsmch * temp end if end if x(j) = ( qtb(j) - sum2 ) / temp end do ! ! Test whether the Gauss-Newton direction is acceptable. ! wa1(1:n) = 0.0D+00 wa2(1:n) = diag(1:n) * x(1:n) qnorm = enorm ( n, wa2 ) if ( qnorm <= delta ) then return end if ! ! The Gauss-Newton direction is not acceptable. ! Calculate the scaled gradient direction. ! l = 1 do j = 1, n temp = qtb(j) do i = j, n wa1(i) = wa1(i) + r(l) * temp l = l + 1 end do wa1(j) = wa1(j) / diag(j) end do ! ! Calculate the norm of the scaled gradient. ! Test for the special case in which the scaled gradient is zero. ! gnorm = enorm ( n, wa1 ) sgnorm = 0.0D+00 alpha = delta / qnorm if ( gnorm /= 0.0D+00 ) then ! ! Calculate the point along the scaled gradient which minimizes the quadratic. ! wa1(1:n) = ( wa1(1:n) / gnorm ) / diag(1:n) l = 1 do j = 1, n sum2 = 0.0D+00 do i = j, n sum2 = sum2 + r(l) * wa1(i) l = l + 1 end do wa2(j) = sum2 end do temp = enorm ( n, wa2 ) sgnorm = ( gnorm / temp ) / temp ! ! Test whether the scaled gradient direction is acceptable. ! alpha = 0.0D+00 ! ! The scaled gradient direction is not acceptable. ! Calculate the point along the dogleg at which the quadratic is minimized. ! if ( sgnorm < delta ) then bnorm = enorm ( n, qtb ) temp = ( bnorm / gnorm ) * ( bnorm / qnorm ) * ( sgnorm / delta ) temp = temp - ( delta / qnorm ) * ( sgnorm / delta) ** 2 & + sqrt ( ( temp - ( delta / qnorm ) ) ** 2 & + ( 1.0D+00 - ( delta / qnorm ) ** 2 ) & * ( 1.0D+00 - ( sgnorm / delta ) ** 2 ) ) alpha = ( ( delta / qnorm ) * ( 1.0D+00 - ( sgnorm / delta ) ** 2 ) ) & / temp end if end if ! ! Form appropriate convex combination of the Gauss-Newton ! direction and the scaled gradient direction. ! temp = ( 1.0D+00 - alpha ) * min ( sgnorm, delta ) x(1:n) = temp * wa1(1:n) + alpha * x(1:n) return endsubroutine dogleg function enorm( n, x ) !*****************************************************************************80 ! !! ENORM computes the Euclidean norm of a vector. ! ! Discussion: ! ! This is an extremely simplified version of the original ENORM ! routine, which has been renamed to "ENORM2". ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, is the length of the vector. ! ! Input, real ( kind = 8 ) X(N), the vector whose norm is desired. ! ! Output, real ( kind = 8 ) ENORM, the Euclidean norm of the vector. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) x(n) real ( kind = 8 ) enorm enorm = sqrt ( sum ( x(1:n) ** 2 )) return endfunction enorm function enorm2 ( n, x ) !*****************************************************************************80 ! !! ENORM2 computes the Euclidean norm of a vector. ! ! Discussion: ! ! This routine was named ENORM. It has been renamed "ENORM2", ! and a simplified routine has been substituted. ! ! The Euclidean norm is computed by accumulating the sum of ! squares in three different sums. The sums of squares for the ! small and large components are scaled so that no overflows ! occur. Non-destructive underflows are permitted. Underflows ! and overflows do not occur in the computation of the unscaled ! sum of squares for the intermediate components. ! ! The definitions of small, intermediate and large components ! depend on two constants, RDWARF and RGIANT. The main ! restrictions on these constants are that RDWARF^2 not ! underflow and RGIANT^2 not overflow. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1 ! Argonne National Laboratory, ! Argonne, Illinois. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, is the length of the vector. ! ! Input, real ( kind = 8 ) X(N), the vector whose norm is desired. ! ! Output, real ( kind = 8 ) ENORM2, the Euclidean norm of the vector. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) agiant real ( kind = 8 ) enorm2 integer ( kind = 4 ) i real ( kind = 8 ) rdwarf real ( kind = 8 ) rgiant real ( kind = 8 ) s1 real ( kind = 8 ) s2 real ( kind = 8 ) s3 real ( kind = 8 ) x(n) real ( kind = 8 ) xabs real ( kind = 8 ) x1max real ( kind = 8 ) x3max rdwarf = sqrt ( tiny ( rdwarf ) ) rgiant = sqrt ( huge ( rgiant ) ) s1 = 0.0D+00 s2 = 0.0D+00 s3 = 0.0D+00 x1max = 0.0D+00 x3max = 0.0D+00 agiant = rgiant / real ( n, kind = 8 ) do i = 1, n xabs = abs ( x(i) ) if ( xabs <= rdwarf ) then if ( x3max < xabs ) then s3 = 1.0D+00 + s3 * ( x3max / xabs ) ** 2 x3max = xabs else if ( xabs /= 0.0D+00 ) then s3 = s3 + ( xabs / x3max ) ** 2 end if else if ( agiant <= xabs ) then if ( x1max < xabs ) then s1 = 1.0D+00 + s1 * ( x1max / xabs ) ** 2 x1max = xabs else s1 = s1 + ( xabs / x1max ) ** 2 end if else s2 = s2 + xabs ** 2 end if end do ! ! Calculation of norm. ! if ( s1 /= 0.0D+00 ) then enorm2 = x1max * sqrt ( s1 + ( s2 / x1max ) / x1max ) else if ( s2 /= 0.0D+00 ) then if ( x3max <= s2 ) then enorm2 = sqrt ( s2 * ( 1.0D+00 + ( x3max / s2 ) * ( x3max * s3 ) ) ) else enorm2 = sqrt ( x3max * ( ( s2 / x3max ) + ( x3max * s3 ) ) ) end if else enorm2 = x3max * sqrt ( s3 ) end if return endfunction enorm2 subroutine fdjac1 ( fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn ) !*****************************************************************************80 ! !! FDJAC1 estimates an N by N jacobian matrix using forward differences. ! ! Discussion: ! ! This function computes a forward-difference approximation ! to the N by N jacobian matrix associated with a specified ! problem of N functions in N variables. If the jacobian has ! a banded form, then function evaluations are saved by only ! approximating the nonzero terms. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions. The routine should have the form: ! subroutine fcn ( n, x, fvec, iflag ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) fvec(n) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) N, the number of functions and variables. ! ! Input, real ( kind = 8 ) X(N), the point where the jacobian is evaluated. ! ! Input, real ( kind = 8 ) FVEC(N), the functions evaluated at X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), the N by N approximate ! jacobian matrix. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC, which ! must not be less than N. ! ! Output, integer ( kind = 4 ) IFLAG, is an error flag returned by FCN. ! If FCN returns a nonzero value of IFLAG, then this routine returns ! immediately to the calling program, with the value of IFLAG. ! ! Input, integer ( kind = 4 ) ML, MU, specify the number of subdiagonals and ! superdiagonals within the band of the jacobian matrix. If the ! jacobian is not banded, set ML and MU to N-1. ! ! Input, real ( kind = 8 ) EPSFCN, is 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. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) n real ( kind = 8 ) eps real ( kind = 8 ) epsfcn real ( kind = 8 ) epsmch external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fvec(n) real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) ml integer ( kind = 4 ) msum integer ( kind = 4 ) mu real ( kind = 8 ) temp real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) x(n) epsmch = epsilon ( epsmch ) eps = sqrt ( max ( epsfcn, epsmch ) ) msum = ml + mu + 1 ! ! Computation of dense approximate jacobian. ! if ( n <= msum ) then do j = 1, n temp = x(j) h = eps * abs ( temp ) if ( h == 0.0D+00 ) then h = eps end if iflag = 1 x(j) = temp + h call fcn ( n, x, wa1, iflag ) if ( iflag < 0 ) then exit end if x(j) = temp fjac(1:n,j) = ( wa1(1:n) - fvec(1:n) ) / h end do else ! ! Computation of banded approximate jacobian. ! do k = 1, msum do j = k, n, msum wa2(j) = x(j) h = eps * abs ( wa2(j) ) if ( h == 0.0D+00 ) then h = eps end if x(j) = wa2(j) + h end do iflag = 1 call fcn ( n, x, wa1, iflag ) if ( iflag < 0 ) then exit end if do j = k, n, msum x(j) = wa2(j) h = eps * abs ( wa2(j) ) if ( h == 0.0D+00 ) then h = eps end if fjac(1:n,j) = 0.0D+00 do i = 1, n if ( j - mu <= i .and. i <= j + ml ) then fjac(i,j) = ( wa1(i) - fvec(i) ) / h end if end do end do end do end if return endsubroutine fdjac1 subroutine fdjac2 ( fcn, m, n, x, fvec, fjac, ldfjac, iflag, epsfcn ) !*****************************************************************************80 ! !! FDJAC2 estimates an M by N jacobian matrix using forward differences. ! ! Discussion: ! ! This function computes a forward-difference approximation ! to the M by N jacobian matrix associated with a specified ! problem of M functions in N variables. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions. The routine should have the form: ! subroutine fcn ( m, n, x, fvec, iflag ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) M, is the number of functions. ! ! Input, integer ( kind = 4 ) N, is the number of variables. ! N must not exceed M. ! ! Input, real ( kind = 8 ) X(N), the point where the jacobian is evaluated. ! ! Input, real ( kind = 8 ) FVEC(M), the functions evaluated at X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), the M by N approximate ! jacobian matrix. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC, ! which must not be less than M. ! ! Output, integer ( kind = 4 ) IFLAG, is an error flag returned by FCN. ! If FCN returns a nonzero value of IFLAG, then this routine returns ! immediately to the calling program, with the value of IFLAG. ! ! Input, real ( kind = 8 ) EPSFCN, is 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. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) eps real ( kind = 8 ) epsfcn real ( kind = 8 ) epsmch external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fvec(m) real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) j real ( kind = 8 ) temp real ( kind = 8 ) wa(m) real ( kind = 8 ) x(n) epsmch = epsilon ( epsmch ) eps = sqrt ( max ( epsfcn, epsmch ) ) do j = 1, n temp = x(j) h = eps * abs ( temp ) if ( h == 0.0D+00 ) then h = eps end if iflag = 1 x(j) = temp + h call fcn ( m, n, x, wa, iflag ) if ( iflag < 0 ) then exit end if x(j) = temp fjac(1:m,j) = ( wa(1:m) - fvec(1:m) ) / h end do return endsubroutine fdjac2 subroutine hybrd ( fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, & factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf ) !*****************************************************************************80 ! !! HYBRD seeks a zero of N nonlinear equations in N variables. ! ! Discussion: ! ! HYBRD finds 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. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions. The routine should have the form: ! subroutine fcn ( n, x, fvec, iflag ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) fvec(n) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) N, the number of functions and variables. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X. ! ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error ! between two consecutive iterates is at most XTOL. XTOL should be ! nonnegative. ! ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of ! calls to FCN is at least MAXFEV by the end of an iteration. ! ! Input, integer ( kind = 4 ) ML, MU, specify the number of subdiagonals and ! superdiagonals within the band of the jacobian matrix. If the jacobian ! is not banded, set ML and MU to at least n - 1. ! ! Input, real ( kind = 8 ) EPSFCN, is 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. ! ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set ! internally. If MODE = 2, then DIAG must contain positive entries that ! serve as multiplicative scale factors for the variables. ! ! Input, integer ( kind = 4 ) MODE, scaling option. ! 1, variables will be scaled internally. ! 2, scaling is specified by the input DIAG vector. ! ! Input, real ( kind = 8 ) FACTOR, determines 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 (0.1, 100) with 100 the recommended value. ! ! Input, integer ( kind = 4 ) NPRINT, 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. ! ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated ! execution, INFO is set to the (negative) value of IFLAG. ! See the description of FCN. ! Otherwise, INFO is set as follows: ! 0, improper input parameters. ! 1, relative error between two consecutive iterates is at most XTOL. ! 2, number of calls to FCN has reached or exceeded MAXFEV. ! 3, XTOL is too small. No further improvement in the approximate ! solution X is possible. ! 4, iteration is not making good progress, as measured by the improvement ! from the last five jacobian evaluations. ! 5, iteration is not making good progress, as measured by the improvement ! from the last ten iterations. ! ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array which contains ! the orthogonal matrix Q produced by the QR factorization of the final ! approximate jacobian. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least N. ! ! Output, real ( kind = 8 ) R(LR), the upper triangular matrix produced by ! the QR factorization of the final approximate jacobian, stored rowwise. ! ! Input, integer ( kind = 4 ) LR, the size of the R array, which must be no ! less than (N*(N+1))/2. ! ! Output, real ( kind = 8 ) QTF(N), contains the vector Q'*FVEC. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) lr integer ( kind = 4 ) n real ( kind = 8 ) actred real ( kind = 8 ) delta real ( kind = 8 ) diag(n) !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsfcn real ( kind = 8 ) epsmch real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fnorm real ( kind = 8 ) fnorm1 real ( kind = 8 ) fvec(n) integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) info integer ( kind = 4 ) iter integer ( kind = 4 ) iwa(1) integer ( kind = 4 ) j logical jeval integer ( kind = 4 ) l integer ( kind = 4 ) maxfev integer ( kind = 4 ) ml integer ( kind = 4 ) mode integer ( kind = 4 ) msum integer ( kind = 4 ) mu integer ( kind = 4 ) ncfail integer ( kind = 4 ) nslow1 integer ( kind = 4 ) nslow2 integer ( kind = 4 ) ncsuc integer ( kind = 4 ) nfev integer ( kind = 4 ) nprint logical pivot real ( kind = 8 ) pnorm real ( kind = 8 ) prered real ( kind = 8 ) qtf(n) real ( kind = 8 ) r(lr) real ( kind = 8 ) ratio logical sing real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) wa3(n) real ( kind = 8 ) wa4(n) real ( kind = 8 ) x(n) real ( kind = 8 ) xnorm real ( kind = 8 ) xtol epsmch = epsilon ( epsmch ) info = 0 iflag = 0 nfev = 0 ! ! Check the input parameters for errors. ! if ( n <= 0 ) then return else if ( xtol < 0.0D+00 ) then return else if ( maxfev <= 0 ) then return else if ( ml < 0 ) then return else if ( mu < 0 ) then return else if ( factor <= 0.0D+00 ) then return else if ( ldfjac < n ) then return else if ( lr < ( n * ( n + 1 ) ) / 2 ) then return end if if ( mode == 2 ) then do j = 1, n if ( diag(j) <= 0.0D+00 ) then go to 300 end if end do end if ! ! Evaluate the function at the starting point ! and calculate its norm. ! iflag = 1 call fcn ( n, x, fvec, iflag ) nfev = 1 if ( iflag < 0 ) then go to 300 end if fnorm = enorm ( n, fvec ) ! ! Determine the number of calls to FCN needed to compute the jacobian matrix. ! msum = min ( ml + mu + 1, n ) ! ! Initialize iteration counter and monitors. ! iter = 1 ncsuc = 0 ncfail = 0 nslow1 = 0 nslow2 = 0 ! ! Beginning of the outer loop. ! 30 continue jeval = .true. ! ! Calculate the jacobian matrix. ! iflag = 2 call fdjac1 ( fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn ) nfev = nfev + msum if ( iflag < 0 ) then go to 300 end if ! ! Compute the QR factorization of the jacobian. ! pivot = .false. call qrfac ( n, n, fjac, ldfjac, pivot, iwa, 1, wa1, wa2 ) ! ! On the first iteration, if MODE is 1, scale according ! to the norms of the columns of the initial jacobian. ! if ( iter == 1 ) then if ( mode /= 2 ) then diag(1:n) = wa2(1:n) do j = 1, n if ( wa2(j) == 0.0D+00 ) then diag(j) = 1.0D+00 end if end do end if ! ! On the first iteration, calculate the norm of the scaled X ! and initialize the step bound DELTA. ! wa3(1:n) = diag(1:n) * x(1:n) xnorm = enorm ( n, wa3 ) delta = factor * xnorm if ( delta == 0.0D+00 ) then delta = factor end if end if ! ! Form Q' * FVEC and store in QTF. ! qtf(1:n) = fvec(1:n) do j = 1, n if ( fjac(j,j) /= 0.0D+00 ) then temp = - dot_product ( qtf(j:n), fjac(j:n,j) ) / fjac(j,j) qtf(j:n) = qtf(j:n) + fjac(j:n,j) * temp end if end do ! ! Copy the triangular factor of the QR factorization into R. ! sing = .false. do j = 1, n l = j do i = 1, j - 1 r(l) = fjac(i,j) l = l + n - i end do r(l) = wa1(j) if ( wa1(j) == 0.0D+00 ) then sing = .true. end if end do ! ! Accumulate the orthogonal factor in FJAC. ! call qform ( n, n, fjac, ldfjac ) ! ! Rescale if necessary. ! if ( mode /= 2 ) then do j = 1, n diag(j) = max ( diag(j), wa2(j) ) end do end if ! ! Beginning of the inner loop. ! 180 continue ! ! If requested, call FCN to enable printing of iterates. ! if ( 0 < nprint ) then iflag = 0 if ( mod ( iter - 1, nprint ) == 0 ) then call fcn ( n, x, fvec, iflag ) end if if ( iflag < 0 ) then go to 300 end if end if ! ! Determine the direction P. ! call dogleg ( n, r, lr, diag, qtf, delta, wa1 ) ! ! Store the direction P and X + P. ! Calculate the norm of P. ! wa1(1:n) = - wa1(1:n) wa2(1:n) = x(1:n) + wa1(1:n) wa3(1:n) = diag(1:n) * wa1(1:n) pnorm = enorm ( n, wa3 ) ! ! On the first iteration, adjust the initial step bound. ! if ( iter == 1 ) then delta = min ( delta, pnorm ) end if ! ! 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 go to 300 end if fnorm1 = enorm ( n, wa4 ) ! ! Compute the scaled actual reduction. ! actred = -1.0D+00 if ( fnorm1 < fnorm ) then actred = 1.0D+00 - ( fnorm1 / fnorm ) ** 2 endif ! ! Compute the scaled predicted reduction. ! l = 1 do i = 1, n sum2 = 0.0D+00 do j = i, n sum2 = sum2 + r(l) * wa1(j) l = l + 1 end do wa3(i) = qtf(i) + sum2 end do temp = enorm ( n, wa3 ) prered = 0.0D+00 if ( temp < fnorm ) then prered = 1.0D+00 - ( temp / fnorm ) ** 2 end if ! ! Compute the ratio of the actual to the predicted reduction. ! ratio = 0.0D+00 if ( 0.0D+00 < prered ) then ratio = actred / prered end if ! ! Update the step bound. ! if ( ratio < 0.1D+00 ) then ncsuc = 0 ncfail = ncfail + 1 delta = 0.5D+00 * delta else ncfail = 0 ncsuc = ncsuc + 1 if ( 0.5D+00 <= ratio .or. 1 < ncsuc ) then delta = max ( delta, pnorm / 0.5D+00 ) end if if ( abs ( ratio - 1.0D+00 ) <= 0.1D+00 ) then delta = pnorm / 0.5D+00 end if end if ! ! Test for successful iteration. ! ! Successful iteration. ! Update X, FVEC, and their norms. ! if ( 0.0001D+00 <= ratio ) then x(1:n) = wa2(1:n) wa2(1:n) = diag(1:n) * x(1:n) fvec(1:n) = wa4(1:n) xnorm = enorm ( n, wa2 ) fnorm = fnorm1 iter = iter + 1 end if ! ! Determine the progress of the iteration. ! nslow1 = nslow1 + 1 if ( 0.001D+00 <= actred ) then nslow1 = 0 end if if ( jeval ) then nslow2 = nslow2 + 1 end if if ( 0.1D+00 <= actred ) then nslow2 = 0 end if ! ! Test for convergence. ! if ( delta <= xtol * xnorm .or. fnorm == 0.0D+00 ) then info = 1 end if if ( info /= 0 ) then go to 300 end if ! ! Tests for termination and stringent tolerances. ! if ( maxfev <= nfev ) then info = 2 end if if ( 0.1D+00 * max ( 0.1D+00 * delta, pnorm ) <= epsmch * xnorm ) then info = 3 end if if ( nslow2 == 5 ) then info = 4 end if if ( nslow1 == 10 ) then info = 5 end if if ( info /= 0 ) then go to 300 end if ! ! Criterion for recalculating jacobian approximation ! by forward differences. ! if ( ncfail == 2 ) then go to 290 end if ! ! Calculate the rank one modification to the jacobian ! and update QTF if necessary. ! do j = 1, n sum2 = dot_product ( wa4(1:n), fjac(1:n,j) ) wa2(j) = ( sum2 - wa3(j) ) / pnorm wa1(j) = diag(j) * ( ( diag(j) * wa1(j) ) / pnorm ) if ( 0.0001D+00 <= ratio ) then qtf(j) = sum2 end if end do ! ! 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. go to 180 290 continue ! ! End of the outer loop. ! go to 30 300 continue ! ! Termination, either normal or user imposed. ! if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, iflag ) end if return endsubroutine hybrd subroutine hybrd1 ( fcn, n, x, fvec, tol, info ) !*****************************************************************************80 ! !! HYBRD1 seeks a zero of N nonlinear equations in N variables. ! ! Discussion: ! ! HYBRD1 finds a zero of a system of N nonlinear functions in N variables ! by a modification of the Powell hybrid method. This is done by using the ! more general nonlinear equation solver HYBRD. The user must provide a ! subroutine which calculates the functions. The jacobian is then ! calculated by a forward-difference approximation. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 19 August 2016 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions. The routine should have the form: ! subroutine fcn ( n, x, fvec, iflag ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) fvec(n) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) N, the number of functions and variables. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X. ! ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm ! estimates that the relative error between X and the solution is at ! most TOL. TOL should be nonnegative. ! ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated ! execution, INFO is set to the (negative) value of IFLAG. See the ! description of FCN. ! Otherwise, INFO is set as follows: ! 0, improper input parameters. ! 1, algorithm estimates that the relative error between X and the ! solution is at most TOL. ! 2, number of calls to FCN has reached or exceeded 200*(N+1). ! 3, TOL is too small. No further improvement in the approximate ! solution X is possible. ! 4, the iteration is not making good progress. ! implicit none integer ( kind = 4 ) lwa integer ( kind = 4 ) n real ( kind = 8 ) diag(n) real ( kind = 8 ) epsfcn real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(n,n) real ( kind = 8 ) fvec(n) integer ( kind = 4 ) info integer ( kind = 4 ) j integer ( kind = 4 ) ldfjac integer ( kind = 4 ) lr integer ( kind = 4 ) maxfev integer ( kind = 4 ) ml integer ( kind = 4 ) mode integer ( kind = 4 ) mu integer ( kind = 4 ) nfev integer ( kind = 4 ) nprint real ( kind = 8 ) qtf(n) real ( kind = 8 ) r((n*(n+1))/2) real ( kind = 8 ) tol real ( kind = 8 ) x(n) real ( kind = 8 ) xtol if ( n <= 0 ) then info = 0 return end if if ( tol < 0.0D+00 ) then info = 0 return end if xtol = tol maxfev = 200 * ( n + 1 ) ml = n - 1 mu = n - 1 epsfcn = 0.0D+00 diag(1:n) = 1.0D+00 mode = 2 factor = 100.0D+00 nprint = 0 info = 0 nfev = 0 fjac(1:n,1:n) = 0.0D+00 ldfjac = n r(1:(n*(n+1))/2) = 0.0D+00 lr = ( n * ( n + 1 ) ) / 2 qtf(1:n) = 0.0D+00 call hybrd ( fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, & factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf ) if ( info == 5 ) then info = 4 end if return endsubroutine hybrd1 subroutine hybrj ( fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, & factor, nprint, info, nfev, njev, r, lr, qtf ) !*****************************************************************************80 ! !! HYBRJ seeks a zero of N nonlinear equations in N variables. ! ! Discussion: ! ! HYBRJ finds 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. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions and the jacobian. FCN should have the form: ! ! subroutine fcn ( n, x, fvec, fjac, ldfjac, iflag ) ! integer ( kind = 4 ) ldfjac ! integer ( kind = 4 ) n ! real ( kind = 8 ) fjac(ldfjac,n) ! real ( kind = 8 ) fvec(n) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and ! return this matrix in FJAC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) N, the number of functions and variables. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N matrix, containing ! the orthogonal matrix Q produced by the QR factorization ! of the final approximate jacobian. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of the ! array FJAC. LDFJAC must be at least N. ! ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error ! between two consecutive iterates is at most XTOL. XTOL should be ! nonnegative. ! ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of ! calls to FCN is at least MAXFEV by the end of an iteration. ! ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set ! internally. If MODE = 2, then DIAG must contain positive entries that ! serve as multiplicative scale factors for the variables. ! ! Input, integer ( kind = 4 ) MODE, scaling option. ! 1, variables will be scaled internally. ! 2, scaling is specified by the input DIAG vector. ! ! Input, real ( kind = 8 ) FACTOR, determines 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 (0.1, 100) with 100 the recommended value. ! ! Input, integer ( kind = 4 ) NPRINT, 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. ! ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated ! execution, INFO is set to the (negative) value of IFLAG. ! See the description of FCN. Otherwise, INFO is set as follows: ! 0, improper input parameters. ! 1, relative error between two consecutive iterates is at most XTOL. ! 2, number of calls to FCN with IFLAG = 1 has reached MAXFEV. ! 3, XTOL is too small. No further improvement in ! the approximate solution X is possible. ! 4, iteration is not making good progress, as measured by the ! improvement from the last five jacobian evaluations. ! 5, iteration is not making good progress, as measured by the ! improvement from the last ten iterations. ! ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN ! with IFLAG = 1. ! ! Output, integer ( kind = 4 ) NJEV, the number of calls to FCN ! with IFLAG = 2. ! ! Output, real ( kind = 8 ) R(LR), the upper triangular matrix produced ! by the QR factorization of the final approximate jacobian, stored rowwise. ! ! Input, integer ( kind = 4 ) LR, the size of the R array, which must ! be no less than (N*(N+1))/2. ! ! Output, real ( kind = 8 ) QTF(N), contains the vector Q'*FVEC. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) lr integer ( kind = 4 ) n real ( kind = 8 ) actred real ( kind = 8 ) delta real ( kind = 8 ) diag(n) !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsmch real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fnorm real ( kind = 8 ) fnorm1 real ( kind = 8 ) fvec(n) integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) info integer ( kind = 4 ) iter integer ( kind = 4 ) iwa(1) integer ( kind = 4 ) j logical jeval integer ( kind = 4 ) l integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) ncfail integer ( kind = 4 ) nslow1 integer ( kind = 4 ) nslow2 integer ( kind = 4 ) ncsuc integer ( kind = 4 ) nfev integer ( kind = 4 ) njev integer ( kind = 4 ) nprint logical pivot real ( kind = 8 ) pnorm real ( kind = 8 ) prered real ( kind = 8 ) qtf(n) real ( kind = 8 ) r(lr) real ( kind = 8 ) ratio logical sing real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) wa3(n) real ( kind = 8 ) wa4(n) real ( kind = 8 ) x(n) real ( kind = 8 ) xnorm real ( kind = 8 ) xtol epsmch = epsilon ( epsmch ) info = 0 iflag = 0 nfev = 0 njev = 0 ! ! Check the input parameters for errors. ! if ( n <= 0 ) then if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if if ( ldfjac < n .or. & xtol < 0.0D+00 .or. & maxfev <= 0 .or. & factor <= 0.0D+00 .or. & lr < ( n * ( n + 1 ) ) / 2 ) then if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if if ( mode == 2 ) then do j = 1, n if ( diag(j) <= 0.0D+00 ) then if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if end do end if ! ! 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 ) then if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if fnorm = enorm ( n, fvec ) ! ! Initialize iteration counter and monitors. ! iter = 1 ncsuc = 0 ncfail = 0 nslow1 = 0 nslow2 = 0 ! ! Beginning of the outer loop. ! do jeval = .true. ! ! Calculate the jacobian matrix. ! iflag = 2 call fcn ( n, x, fvec, fjac, ldfjac, iflag ) njev = njev + 1 if ( iflag < 0 ) then info = iflag iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if ! ! Compute the QR factorization of the jacobian. ! pivot = .false. call qrfac ( n, n, fjac, ldfjac, pivot, iwa, 1, wa1, wa2 ) ! ! On the first iteration, if MODE is 1, scale according ! to the norms of the columns of the initial jacobian. ! if ( iter == 1 ) then if ( mode /= 2 ) then diag(1:n) = wa2(1:n) do j = 1, n if ( wa2(j) == 0.0D+00 ) then diag(j) = 1.0D+00 end if end do end if ! ! On the first iteration, calculate the norm of the scaled X ! and initialize the step bound DELTA. ! wa3(1:n) = diag(1:n) * x(1:n) xnorm = enorm ( n, wa3 ) delta = factor * xnorm if ( delta == 0.0D+00 ) then delta = factor end if end if ! ! Form Q'*FVEC and store in QTF. ! qtf(1:n) = fvec(1:n) do j = 1, n if ( fjac(j,j) /= 0.0D+00 ) then sum2 = 0.0D+00 do i = j, n sum2 = sum2 + fjac(i,j) * qtf(i) end do temp = - sum2 / fjac(j,j) do i = j, n qtf(i) = qtf(i) + fjac(i,j) * temp end do end if end do ! ! Copy the triangular factor of the QR factorization into R. ! sing = .false. do j = 1, n l = j do i = 1, j - 1 r(l) = fjac(i,j) l = l + n - i end do r(l) = wa1(j) if ( wa1(j) == 0.0D+00 ) then sing = .true. end if end do ! ! Accumulate the orthogonal factor in FJAC. ! call qform ( n, n, fjac, ldfjac ) ! ! Rescale if necessary. ! if ( mode /= 2 ) then do j = 1, n diag(j) = max ( diag(j), wa2(j) ) end do end if ! ! Beginning of the inner loop. ! do ! ! If requested, call FCN to enable printing of iterates. ! if ( 0 < nprint ) then iflag = 0 if ( mod ( iter - 1, nprint ) == 0 ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if if ( iflag < 0 ) then info = iflag iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if end if ! ! Determine the direction P. ! call dogleg ( n, r, lr, diag, qtf, delta, wa1 ) ! ! Store the direction P and X + P. ! Calculate the norm of P. ! wa1(1:n) = - wa1(1:n) wa2(1:n) = x(1:n) + wa1(1:n) wa3(1:n) = diag(1:n) * wa1(1:n) pnorm = enorm ( n, wa3 ) ! ! On the first iteration, adjust the initial step bound. ! if ( iter == 1 ) then delta = min ( delta, pnorm ) end if ! ! 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 info = iflag iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if fnorm1 = enorm ( n, wa4 ) ! ! Compute the scaled actual reduction. ! actred = -1.0D+00 if ( fnorm1 < fnorm ) then actred = 1.0D+00 - ( fnorm1 / fnorm ) ** 2 end if ! ! Compute the scaled predicted reduction. ! l = 1 do i = 1, n sum2 = 0.0D+00 do j = i, n sum2 = sum2 + r(l) * wa1(j) l = l + 1 end do wa3(i) = qtf(i) + sum2 end do temp = enorm ( n, wa3 ) prered = 0.0D+00 if ( temp < fnorm ) then prered = 1.0D+00 - ( temp / fnorm ) ** 2 end if ! ! Compute the ratio of the actual to the predicted reduction. ! if ( 0.0D+00 < prered ) then ratio = actred / prered else ratio = 0.0D+00 end if ! ! Update the step bound. ! if ( ratio < 0.1D+00 ) then ncsuc = 0 ncfail = ncfail + 1 delta = 0.5D+00 * delta else ncfail = 0 ncsuc = ncsuc + 1 if ( 0.5D+00 <= ratio .or. 1 < ncsuc ) then delta = max ( delta, pnorm / 0.5D+00 ) end if if ( abs ( ratio - 1.0D+00 ) <= 0.1D+00 ) then delta = pnorm / 0.5D+00 end if end if ! ! Test for successful iteration. ! ! ! Successful iteration. ! Update X, FVEC, and their norms. ! if ( 0.0001D+00 <= ratio ) then x(1:n) = wa2(1:n) wa2(1:n) = diag(1:n) * x(1:n) fvec(1:n) = wa4(1:n) xnorm = enorm ( n, wa2 ) fnorm = fnorm1 iter = iter + 1 end if ! ! Determine the progress of the iteration. ! nslow1 = nslow1 + 1 if ( 0.001D+00 <= actred ) then nslow1 = 0 end if if ( jeval ) then nslow2 = nslow2 + 1 end if if ( 0.1D+00 <= actred ) then nslow2 = 0 end if ! ! Test for convergence. ! if ( delta <= xtol * xnorm .or. fnorm == 0.0D+00 ) then info = 1 end if if ( info /= 0 ) then iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if ! ! Tests for termination and stringent tolerances. ! if ( maxfev <= nfev ) then info = 2 end if if ( 0.1D+00 * max ( 0.1D+00 * delta, pnorm ) <= epsmch * xnorm ) then info = 3 end if if ( nslow2 == 5 ) then info = 4 end if if ( nslow1 == 10 ) then info = 5 end if if ( info /= 0 ) then iflag = 0 if ( 0 < nprint ) then call fcn ( n, x, fvec, fjac, ldfjac, iflag ) end if return end if ! ! Criterion for recalculating jacobian. ! if ( ncfail == 2 ) then exit end if ! ! Calculate the rank one modification to the jacobian ! and update QTF if necessary. ! do j = 1, n sum2 = dot_product ( wa4(1:n), fjac(1:n,j) ) wa2(j) = ( sum2 - wa3(j) ) / pnorm wa1(j) = diag(j) * ( ( diag(j) * wa1(j) ) / pnorm ) if ( 0.0001D+00 <= ratio ) then qtf(j) = sum2 end if end do ! ! 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 do ! ! End of the outer loop. ! end do endsubroutine hybrj subroutine hybrj1 ( fcn, n, x, fvec, fjac, ldfjac, tol, info ) !*****************************************************************************80 ! !! HYBRJ1 seeks a zero of N equations in N variables by Powell's method. ! ! Discussion: ! ! HYBRJ1 finds a zero of a system of N nonlinear functions in N variables ! by a modification of the Powell hybrid method. This is done by using the ! more general nonlinear equation solver HYBRJ. The user ! must provide a subroutine which calculates the functions ! and the jacobian. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions and the jacobian. FCN should have the form: ! subroutine fcn ( n, x, fvec, fjac, ldfjac, iflag ) ! integer ( kind = 4 ) ldfjac ! integer ( kind = 4 ) n ! real ( kind = 8 ) fjac(ldfjac,n) ! real ( kind = 8 ) fvec(n) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and ! return this matrix in FJAC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) N, the number of functions and variables. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array which contains ! the orthogonal matrix Q produced by the QR factorization of the final ! approximate jacobian. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least N. ! ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm ! estimates that the relative error between X and the solution is at most ! TOL. TOL should be nonnegative. ! ! Output, integer ( kind = 4 ) INFO, error flag. 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: ! 0, improper input parameters. ! 1, algorithm estimates that the relative error between X and the ! solution is at most TOL. ! 2, number of calls to FCN with IFLAG = 1 has reached 100*(N+1). ! 3, TOL is too small. No further improvement in the approximate ! solution X is possible. ! 4, iteration is not making good progress. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) n real ( kind = 8 ) diag(n) real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fvec(n) integer ( kind = 4 ) info integer ( kind = 4 ) j integer ( kind = 4 ) lr integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) njev integer ( kind = 4 ) nprint real ( kind = 8 ) qtf(n) real ( kind = 8 ) r((n*(n+1))/2) real ( kind = 8 ) tol real ( kind = 8 ) x(n) real ( kind = 8 ) xtol info = 0 if ( n <= 0 ) then return else if ( ldfjac < n ) then return else if ( tol < 0.0D+00 ) then return end if maxfev = 100 * ( n + 1 ) xtol = tol mode = 2 diag(1:n) = 1.0D+00 factor = 100.0D+00 nprint = 0 lr = ( n * ( n + 1 ) ) / 2 call hybrj ( fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, & factor, nprint, info, nfev, njev, r, lr, qtf ) if ( info == 5 ) then info = 4 end if return endsubroutine hybrj1 subroutine lmder ( fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf ) !*****************************************************************************80 ! !! LMDER minimizes M functions in N variables by the Levenberg-Marquardt method. ! ! Discussion: ! ! LMDER minimizes the sum of the squares of M nonlinear functions in ! N variables by a modification of the Levenberg-Marquardt algorithm. ! The user must provide a subroutine which calculates the functions ! and the jacobian. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions and the jacobian. FCN should have the form: ! subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag ) ! integer ( kind = 4 ) ldfjac ! integer ( kind = 4 ) n ! real ( kind = 8 ) fjac(ldfjac,n) ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and ! return this matrix in FJAC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) M, is the number of functions. ! ! Input, integer ( kind = 4 ) N, is the number of variables. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. The upper ! N by N submatrix of FJAC contains an upper triangular matrix R with ! diagonal elements of nonincreasing magnitude such that ! P' * ( JAC' * JAC ) * P = R' * R, ! where P is a permutation matrix and JAC is the final calculated jacobian. ! Column J of P is column IPVT(J) of the identity matrix. The lower ! trapezoidal part of FJAC contains information generated during ! the computation of R. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least M. ! ! Input, real ( kind = 8 ) FTOL. Termination occurs when both the actual ! and predicted relative reductions in the sum of squares are at most FTOL. ! Therefore, FTOL measures the relative error desired in the sum of ! squares. FTOL should be nonnegative. ! ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error ! between two consecutive iterates is at most XTOL. XTOL should be ! nonnegative. ! ! Input, real ( kind = 8 ) GTOL. Termination occurs when the cosine of the ! angle between FVEC and any column of the jacobian is at most GTOL in ! absolute value. Therefore, GTOL measures the orthogonality desired ! between the function vector and the columns of the jacobian. GTOL should ! be nonnegative. ! ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of ! calls to FCN with IFLAG = 1 is at least MAXFEV by the end of an iteration. ! ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set ! internally. If MODE = 2, then DIAG must contain positive entries that ! serve as multiplicative scale factors for the variables. ! ! Input, integer ( kind = 4 ) MODE, scaling option. ! 1, variables will be scaled internally. ! 2, scaling is specified by the input DIAG vector. ! ! Input, real ( kind = 8 ) FACTOR, determines 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 (0.1, 100) with 100 the recommended value. ! ! Input, integer ( kind = 4 ) NPRINT, 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. ! ! Output, integer ( kind = 4 ) INFO, error flag. 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: ! 0, improper input parameters. ! 1, both actual and predicted relative reductions in the sum of ! squares are at most FTOL. ! 2, relative error between two consecutive iterates is at most XTOL. ! 3, conditions for INFO = 1 and INFO = 2 both hold. ! 4, the cosine of the angle between FVEC and any column of the jacobian ! is at most GTOL in absolute value. ! 5, number of calls to FCN with IFLAG = 1 has reached MAXFEV. ! 6, FTOL is too small. No further reduction in the sum of squares ! is possible. ! 7, XTOL is too small. No further improvement in the approximate ! solution X is possible. ! 8, GTOL is too small. FVEC is orthogonal to the columns of the ! jacobian to machine precision. ! ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN with ! IFLAG = 1. ! ! Output, integer ( kind = 4 ) NJEV, the number of calls to FCN with ! IFLAG = 2. ! ! Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P ! such that JAC*P = Q*R, where JAC is the final calculated jacobian, Q is ! orthogonal (not stored), and R is upper triangular with diagonal ! elements of nonincreasing magnitude. Column J of P is column ! IPVT(J) of the identity matrix. ! ! Output, real ( kind = 8 ) QTF(N), contains the first N elements of Q'*FVEC. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) actred real ( kind = 8 ) delta real ( kind = 8 ) diag(n) real ( kind = 8 ) dirder !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsmch real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fnorm real ( kind = 8 ) fnorm1 real ( kind = 8 ) ftol real ( kind = 8 ) fvec(m) real ( kind = 8 ) gnorm real ( kind = 8 ) gtol integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) info integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) iter integer ( kind = 4 ) j integer ( kind = 4 ) l integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) njev integer ( kind = 4 ) nprint real ( kind = 8 ) par logical pivot real ( kind = 8 ) pnorm real ( kind = 8 ) prered real ( kind = 8 ) qtf(n) real ( kind = 8 ) ratio real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) temp1 real ( kind = 8 ) temp2 real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) wa3(n) real ( kind = 8 ) wa4(m) real ( kind = 8 ) xnorm real ( kind = 8 ) x(n) real ( kind = 8 ) xtol epsmch = epsilon ( epsmch ) info = 0 iflag = 0 nfev = 0 njev = 0 ! ! Check the input parameters for errors. ! if ( n <= 0 ) then go to 300 end if if ( m < n ) then go to 300 end if if ( ldfjac < m & .or. ftol < 0.0D+00 .or. xtol < 0.0D+00 .or. gtol < 0.0D+00 & .or. maxfev <= 0 .or. factor <= 0.0D+00 ) then go to 300 end if if ( mode == 2 ) then do j = 1, n if ( diag(j) <= 0.0D+00 ) then go to 300 end if end do end if ! ! Evaluate the function at the starting point and calculate its norm. ! iflag = 1 call fcn ( m, n, x, fvec, fjac, ldfjac, iflag ) nfev = 1 if ( iflag < 0 ) then go to 300 end if fnorm = enorm ( m, fvec ) ! ! Initialize Levenberg-Marquardt parameter and iteration counter. ! par = 0.0D+00 iter = 1 ! ! Beginning of the outer loop. ! 30 continue ! ! Calculate the jacobian matrix. ! iflag = 2 call fcn ( m, n, x, fvec, fjac, ldfjac, iflag ) njev = njev + 1 if ( iflag < 0 ) then go to 300 end if ! ! If requested, call FCN to enable printing of iterates. ! if ( 0 < nprint ) then iflag = 0 if ( mod ( iter - 1, nprint ) == 0 ) then call fcn ( m, n, x, fvec, fjac, ldfjac, iflag ) end if if ( iflag < 0 ) then go to 300 end if end if ! ! Compute the QR factorization of the jacobian. ! pivot = .true. call qrfac ( m, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2 ) ! ! 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 diag(1:n) = wa2(1:n) do j = 1, n if ( wa2(j) == 0.0D+00 ) then diag(j) = 1.0D+00 end if end do end if ! ! On the first iteration, calculate the norm of the scaled X ! and initialize the step bound DELTA. ! wa3(1:n) = diag(1:n) * x(1:n) xnorm = enorm ( n, wa3 ) delta = factor * xnorm if ( delta == 0.0D+00 ) then delta = factor end if end if ! ! Form Q'*FVEC and store the first N components in QTF. ! wa4(1:m) = fvec(1:m) do j = 1, n if ( fjac(j,j) /= 0.0D+00 ) then sum2 = dot_product ( wa4(j:m), fjac(j:m,j) ) temp = - sum2 / fjac(j,j) wa4(j:m) = wa4(j:m) + fjac(j:m,j) * temp end if fjac(j,j) = wa1(j) qtf(j) = wa4(j) end do ! ! Compute the norm of the scaled gradient. ! gnorm = 0.0D+00 if ( fnorm /= 0.0D+00 ) then do j = 1, n l = ipvt(j) if ( wa2(l) /= 0.0D+00 ) then sum2 = dot_product ( qtf(1:j), fjac(1:j,j) ) / fnorm gnorm = max ( gnorm, abs ( sum2 / wa2(l) ) ) end if end do end if ! ! Test for convergence of the gradient norm. ! if ( gnorm <= gtol ) then info = 4 go to 300 end if ! ! Rescale if necessary. ! if ( mode /= 2 ) then do j = 1, n diag(j) = max ( diag(j), wa2(j) ) end do end if ! ! Beginning of the inner loop. ! 200 continue ! ! Determine the Levenberg-Marquardt parameter. ! call lmpar ( n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2 ) ! ! Store the direction p and x + p. calculate the norm of p. ! wa1(1:n) = - wa1(1:n) wa2(1:n) = x(1:n) + wa1(1:n) wa3(1:n) = diag(1:n) * wa1(1:n) pnorm = enorm ( n, wa3 ) ! ! On the first iteration, adjust the initial step bound. ! if ( iter == 1 ) then delta = min ( delta, pnorm ) end if ! ! Evaluate the function at x + p and calculate its norm. ! iflag = 1 call fcn ( m, n, wa2, wa4, fjac, ldfjac, iflag ) nfev = nfev + 1 if ( iflag < 0 ) then go to 300 end if fnorm1 = enorm ( m, wa4 ) ! ! Compute the scaled actual reduction. ! actred = -1.0D+00 if ( 0.1D+00 * fnorm1 < fnorm ) then actred = 1.0D+00 - ( fnorm1 / fnorm ) ** 2 end if ! ! Compute the scaled predicted reduction and ! the scaled directional derivative. ! do j = 1, n wa3(j) = 0.0D+00 l = ipvt(j) temp = wa1(l) wa3(1:j) = wa3(1:j) + fjac(1:j,j) * temp end do temp1 = enorm ( n, wa3 ) / fnorm temp2 = ( sqrt ( par ) * pnorm ) / fnorm prered = temp1 ** 2 + temp2 ** 2 / 0.5D+00 dirder = - ( temp1 ** 2 + temp2 ** 2 ) ! ! Compute the ratio of the actual to the predicted reduction. ! if ( prered /= 0.0D+00 ) then ratio = actred / prered else ratio = 0.0D+00 end if ! ! Update the step bound. ! if ( ratio <= 0.25D+00 ) then if ( 0.0D+00 <= actred ) then temp = 0.5D+00 end if if ( actred < 0.0D+00 ) then temp = 0.5D+00 * dirder / ( dirder + 0.5D+00 * actred ) end if if ( 0.1D+00 * fnorm1 >= fnorm .or. temp < 0.1D+00 ) then temp = 0.1D+00 end if delta = temp * min ( delta, pnorm / 0.1D+00 ) par = par / temp else if ( par == 0.0D+00 .or. ratio >= 0.75D+00 ) then delta = 2.0D+00 * pnorm par = 0.5D+00 * par end if end if ! ! Successful iteration. ! ! Update X, FVEC, and their norms. ! if ( 0.0001D+00 <= ratio ) then x(1:n) = wa2(1:n) wa2(1:n) = diag(1:n) * x(1:n) fvec(1:m) = wa4(1:m) xnorm = enorm ( n, wa2 ) fnorm = fnorm1 iter = iter + 1 end if ! ! Tests for convergence. ! if ( abs ( actred) <= ftol .and. & prered <= ftol .and. & 0.5D+00 * ratio <= 1.0D+00 ) then info = 1 end if if ( delta <= xtol * xnorm ) then info = 2 end if if ( abs ( actred) <= ftol .and. prered <= ftol & .and. 0.5D+00 * ratio <= 1.0D+00 .and. info == 2 ) then info = 3 end if if ( info /= 0 ) then go to 300 end if ! ! Tests for termination and stringent tolerances. ! if ( nfev >= maxfev ) then info = 5 end if if ( abs ( actred ) <= epsmch .and. prered <= epsmch & .and. 0.5D+00 * ratio <= 1.0D+00 ) then info = 6 end if if ( delta <= epsmch * xnorm ) then info = 7 end if if ( gnorm <= epsmch ) then info = 8 end if if ( info /= 0 ) then go to 300 end if ! ! End of the inner loop. repeat if iteration unsuccessful. ! if ( ratio < 0.0001D+00 ) then go to 200 end if ! ! End of the outer loop. ! go to 30 300 continue ! ! Termination, either normal or user imposed. ! if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( m, n, x, fvec, fjac, ldfjac, iflag ) end if return endsubroutine lmder subroutine lmder1 ( fcn, m, n, x, fvec, fjac, ldfjac, tol, info ) !*****************************************************************************80 ! !! LMDER1 minimizes M functions in N variables by Levenberg-Marquardt method. ! ! Discussion: ! ! LMDER1 minimizes the sum of the squares of M nonlinear functions in ! N variables by a modification of the Levenberg-Marquardt algorithm. ! This is done by using the more general least-squares solver LMDER. ! The user must provide a subroutine which calculates the functions ! and the jacobian. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions and the jacobian. FCN should have the form: ! subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag ) ! integer ( kind = 4 ) ldfjac ! integer ( kind = 4 ) n ! real ( kind = 8 ) fjac(ldfjac,n) ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and ! return this matrix in FJAC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) M, the number of functions. ! ! Input, integer ( kind = 4 ) N, is the number of variables. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. The upper ! N by N submatrix contains an upper triangular matrix R with ! diagonal elements of nonincreasing magnitude such that ! P' * ( JAC' * JAC ) * P = R' * R, ! where P is a permutation matrix and JAC is the final calculated ! jacobian. Column J of P is column IPVT(J) of the identity matrix. ! The lower trapezoidal part of FJAC contains information generated during ! the computation of R. ! ! Input, integer ( kind = 4 ) LDFJAC, is the leading dimension of FJAC, ! which must be no less than M. ! ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm ! estimates either that the relative error in the sum of squares is at ! most TOL or that the relative error between X and the solution is at ! most TOL. ! ! Output, integer ( kind = 4 ) INFO, error flag. 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: ! 0, improper input parameters. ! 1, algorithm estimates that the relative error in the sum of squares ! is at most TOL. ! 2, algorithm estimates that the relative error between X and the ! solution is at most TOL. ! 3, conditions for INFO = 1 and INFO = 2 both hold. ! 4, FVEC is orthogonal to the columns of the jacobian to machine precision. ! 5, number of calls to FCN with IFLAG = 1 has reached 100*(N+1). ! 6, TOL is too small. No further reduction in the sum of squares is ! possible. ! 7, TOL is too small. No further improvement in the approximate ! solution X is possible. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) diag(n) real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) ftol real ( kind = 8 ) fvec(m) real ( kind = 8 ) gtol integer ( kind = 4 ) info integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) njev integer ( kind = 4 ) nprint real ( kind = 8 ) qtf(n) real ( kind = 8 ) tol real ( kind = 8 ) x(n) real ( kind = 8 ) xtol info = 0 if ( n <= 0 ) then return else if ( m < n ) then return else if ( ldfjac < m ) then return else if ( tol < 0.0D+00 ) then return end if factor = 100.0D+00 maxfev = 100 * ( n + 1 ) ftol = tol xtol = tol gtol = 0.0D+00 mode = 1 nprint = 0 call lmder ( fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf ) if ( info == 8 ) then info = 4 end if return endsubroutine lmder1 subroutine lmdif ( fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, & diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf ) !*****************************************************************************80 ! !! LMDIF minimizes M functions in N variables by the Levenberg-Marquardt method. ! ! Discussion: ! ! LMDIF minimizes the sum of the squares of M nonlinear functions in ! N variables by a modification of the Levenberg-Marquardt algorithm. ! The user must provide a subroutine which calculates the functions. ! The jacobian is then calculated by a forward-difference approximation. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions. The routine should have the form: ! subroutine fcn ( m, n, x, fvec, iflag ) ! integer ( kind = 4 ) m ! integer ( kind = 4 ) n ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) M, the number of functions. ! ! Input, integer ( kind = 4 ) N, the number of variables. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X. ! ! Input, real ( kind = 8 ) FTOL. Termination occurs when both the actual ! and predicted relative reductions in the sum of squares are at most FTOL. ! Therefore, FTOL measures the relative error desired in the sum of ! squares. FTOL should be nonnegative. ! ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error ! between two consecutive iterates is at most XTOL. Therefore, XTOL ! measures the relative error desired in the approximate solution. XTOL ! should be nonnegative. ! ! Input, real ( kind = 8 ) GTOL. termination occurs when the cosine of the ! angle between FVEC and any column of the jacobian is at most GTOL in ! absolute value. Therefore, GTOL measures the orthogonality desired ! between the function vector and the columns of the jacobian. GTOL should ! be nonnegative. ! ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of ! calls to FCN is at least MAXFEV by the end of an iteration. ! ! Input, real ( kind = 8 ) EPSFCN, is 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. ! ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set ! internally. If MODE = 2, then DIAG must contain positive entries that ! serve as multiplicative scale factors for the variables. ! ! Input, integer ( kind = 4 ) MODE, scaling option. ! 1, variables will be scaled internally. ! 2, scaling is specified by the input DIAG vector. ! ! Input, real ( kind = 8 ) FACTOR, determines 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 (0.1, 100) with 100 the recommended value. ! ! Input, integer ( kind = 4 ) NPRINT, 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. ! ! Output, integer ( kind = 4 ) INFO, error flag. 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: ! 0, improper input parameters. ! 1, both actual and predicted relative reductions in the sum of squares ! are at most FTOL. ! 2, relative error between two consecutive iterates is at most XTOL. ! 3, conditions for INFO = 1 and INFO = 2 both hold. ! 4, the cosine of the angle between FVEC and any column of the jacobian ! is at most GTOL in absolute value. ! 5, number of calls to FCN has reached or exceeded MAXFEV. ! 6, FTOL is too small. No further reduction in the sum of squares ! is possible. ! 7, XTOL is too small. No further improvement in the approximate ! solution X is possible. ! 8, GTOL is too small. FVEC is orthogonal to the columns of the ! jacobian to machine precision. ! ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. The upper ! N by N submatrix of FJAC contains an upper triangular matrix R with ! diagonal elements of nonincreasing magnitude such that ! ! P' * ( JAC' * JAC ) * P = R' * R, ! ! where P is a permutation matrix and JAC is the final calculated jacobian. ! Column J of P is column IPVT(J) of the identity matrix. The lower ! trapezoidal part of FJAC contains information generated during ! the computation of R. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least M. ! ! Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P such ! that JAC * P = Q * R, where JAC is the final calculated jacobian, Q is ! orthogonal (not stored), and R is upper triangular with diagonal ! elements of nonincreasing magnitude. Column J of P is column IPVT(J) ! of the identity matrix. ! ! Output, real ( kind = 8 ) QTF(N), the first N elements of Q'*FVEC. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) actred real ( kind = 8 ) delta real ( kind = 8 ) diag(n) real ( kind = 8 ) dirder !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsfcn real ( kind = 8 ) epsmch real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fnorm real ( kind = 8 ) fnorm1 real ( kind = 8 ) ftol real ( kind = 8 ) fvec(m) real ( kind = 8 ) gnorm real ( kind = 8 ) gtol integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) iter integer ( kind = 4 ) info integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) j integer ( kind = 4 ) l integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) nprint real ( kind = 8 ) par logical pivot real ( kind = 8 ) pnorm real ( kind = 8 ) prered real ( kind = 8 ) qtf(n) real ( kind = 8 ) ratio real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) temp1 real ( kind = 8 ) temp2 real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) wa3(n) real ( kind = 8 ) wa4(m) real ( kind = 8 ) x(n) real ( kind = 8 ) xnorm real ( kind = 8 ) xtol epsmch = epsilon ( epsmch ) info = 0 iflag = 0 nfev = 0 if ( n <= 0 ) then go to 300 else if ( m < n ) then go to 300 else if ( ldfjac < m ) then go to 300 else if ( ftol < 0.0D+00 ) then go to 300 else if ( xtol < 0.0D+00 ) then go to 300 else if ( gtol < 0.0D+00 ) then go to 300 else if ( maxfev <= 0 ) then go to 300 else if ( factor <= 0.0D+00 ) then go to 300 end if if ( mode == 2 ) then do j = 1, n if ( diag(j) <= 0.0D+00 ) then go to 300 end if end do end if ! ! Evaluate the function at the starting point and calculate its norm. ! iflag = 1 call fcn ( m, n, x, fvec, iflag ) nfev = 1 if ( iflag < 0 ) then go to 300 end if fnorm = enorm ( m, fvec ) ! ! Initialize Levenberg-Marquardt parameter and iteration counter. ! par = 0.0D+00 iter = 1 ! ! Beginning of the outer loop. ! 30 continue ! ! Calculate the jacobian matrix. ! iflag = 2 call fdjac2 ( fcn, m, n, x, fvec, fjac, ldfjac, iflag, epsfcn ) nfev = nfev + n if ( iflag < 0 ) then go to 300 end if ! ! If requested, call FCN to enable printing of iterates. ! if ( 0 < nprint ) then iflag = 0 if ( mod ( iter - 1, nprint ) == 0 ) then call fcn ( m, n, x, fvec, iflag ) end if if ( iflag < 0 ) then go to 300 end if end if ! ! Compute the QR factorization of the jacobian. ! pivot = .true. call qrfac ( m, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2 ) ! ! 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 diag(1:n) = wa2(1:n) do j = 1, n if ( wa2(j) == 0.0D+00 ) then diag(j) = 1.0D+00 end if end do end if ! ! On the first iteration, calculate the norm of the scaled X ! and initialize the step bound DELTA. ! wa3(1:n) = diag(1:n) * x(1:n) xnorm = enorm ( n, wa3 ) delta = factor * xnorm if ( delta == 0.0D+00 ) then delta = factor end if end if ! ! Form Q' * FVEC and store the first N components in QTF. ! wa4(1:m) = fvec(1:m) do j = 1, n if ( fjac(j,j) /= 0.0D+00 ) then sum2 = dot_product ( wa4(j:m), fjac(j:m,j) ) temp = - sum2 / fjac(j,j) wa4(j:m) = wa4(j:m) + fjac(j:m,j) * temp end if fjac(j,j) = wa1(j) qtf(j) = wa4(j) end do ! ! Compute the norm of the scaled gradient. ! gnorm = 0.0D+00 if ( fnorm /= 0.0D+00 ) then do j = 1, n l = ipvt(j) if ( wa2(l) /= 0.0D+00 ) then sum2 = 0.0D+00 do i = 1, j sum2 = sum2 + fjac(i,j) * ( qtf(i) / fnorm ) end do gnorm = max ( gnorm, abs ( sum2 / wa2(l) ) ) end if end do end if ! ! Test for convergence of the gradient norm. ! if ( gnorm <= gtol ) then info = 4 go to 300 end if ! ! Rescale if necessary. ! if ( mode /= 2 ) then do j = 1, n diag(j) = max ( diag(j), wa2(j) ) end do end if ! ! Beginning of the inner loop. ! 200 continue ! ! Determine the Levenberg-Marquardt parameter. ! call lmpar ( n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2 ) ! ! Store the direction P and X + P. ! Calculate the norm of P. ! wa1(1:n) = -wa1(1:n) wa2(1:n) = x(1:n) + wa1(1:n) wa3(1:n) = diag(1:n) * wa1(1:n) pnorm = enorm ( n, wa3 ) ! ! On the first iteration, adjust the initial step bound. ! if ( iter == 1 ) then delta = min ( delta, pnorm ) end if ! ! Evaluate the function at X + P and calculate its norm. ! iflag = 1 call fcn ( m, n, wa2, wa4, iflag ) nfev = nfev + 1 if ( iflag < 0 ) then go to 300 end if fnorm1 = enorm ( m, wa4 ) ! ! Compute the scaled actual reduction. ! if ( 0.1D+00 * fnorm1 < fnorm ) then actred = 1.0D+00 - ( fnorm1 / fnorm ) ** 2 else actred = -1.0D+00 end if ! ! Compute the scaled predicted reduction and the scaled directional derivative. ! do j = 1, n wa3(j) = 0.0D+00 l = ipvt(j) temp = wa1(l) wa3(1:j) = wa3(1:j) + fjac(1:j,j) * temp end do temp1 = enorm ( n, wa3 ) / fnorm temp2 = ( sqrt ( par ) * pnorm ) / fnorm prered = temp1 ** 2 + temp2 ** 2 / 0.5D+00 dirder = - ( temp1 ** 2 + temp2 ** 2 ) ! ! Compute the ratio of the actual to the predicted reduction. ! ratio = 0.0D+00 if ( prered /= 0.0D+00 ) then ratio = actred / prered end if ! ! Update the step bound. ! if ( ratio <= 0.25D+00 ) then if ( actred >= 0.0D+00 ) then temp = 0.5D+00 endif if ( actred < 0.0D+00 ) then temp = 0.5D+00 * dirder / ( dirder + 0.5D+00 * actred ) end if if ( 0.1D+00 * fnorm1 >= fnorm .or. temp < 0.1D+00 ) then temp = 0.1D+00 end if delta = temp * min ( delta, pnorm / 0.1D+00 ) par = par / temp else if ( par == 0.0D+00 .or. ratio >= 0.75D+00 ) then delta = 2.0D+00 * pnorm par = 0.5D+00 * par end if end if ! ! Test for successful iteration. ! ! ! Successful iteration. update X, FVEC, and their norms. ! if ( 0.0001D+00 <= ratio ) then x(1:n) = wa2(1:n) wa2(1:n) = diag(1:n) * x(1:n) fvec(1:m) = wa4(1:m) xnorm = enorm ( n, wa2 ) fnorm = fnorm1 iter = iter + 1 end if ! ! Tests for convergence. ! if ( abs ( actred) <= ftol .and. prered <= ftol & .and. 0.5D+00 * ratio <= 1.0D+00 ) then info = 1 end if if ( delta <= xtol * xnorm ) then info = 2 end if if ( abs ( actred) <= ftol .and. prered <= ftol & .and. 0.5D+00 * ratio <= 1.0D+00 .and. info == 2 ) info = 3 if ( info /= 0 ) then go to 300 end if ! ! Tests for termination and stringent tolerances. ! if ( maxfev <= nfev ) then info = 5 end if if ( abs ( actred) <= epsmch .and. prered <= epsmch & .and. 0.5D+00 * ratio <= 1.0D+00 ) then info = 6 end if if ( delta <= epsmch * xnorm ) then info = 7 end if if ( gnorm <= epsmch ) then info = 8 end if if ( info /= 0 ) then go to 300 end if ! ! End of the inner loop. Repeat if iteration unsuccessful. ! if ( ratio < 0.0001D+00 ) then go to 200 end if ! ! End of the outer loop. ! go to 30 300 continue ! ! Termination, either normal or user imposed. ! if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( m, n, x, fvec, iflag ) end if return endsubroutine lmdif subroutine lmdif1 ( fcn, m, n, x, fvec, tol, info ) !*****************************************************************************80 ! !! LMDIF1 minimizes M functions in N variables using Levenberg-Marquardt method. ! ! Discussion: ! ! LMDIF1 minimizes the sum of the squares of M nonlinear functions in ! N variables by a modification of the Levenberg-Marquardt algorithm. ! This is done by using the more general least-squares solver LMDIF. ! The user must provide a subroutine which calculates the functions. ! The jacobian is then calculated by a forward-difference approximation. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions. The routine should have the form: ! subroutine fcn ( m, n, x, fvec, iflag ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! To terminate the algorithm, FCN may set IFLAG negative on return. ! ! Input, integer ( kind = 4 ) M, the number of functions. ! ! Input, integer ( kind = 4 ) N, the number of variables. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X. ! ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm ! estimates either that the relative error in the sum of squares is at ! most TOL or that the relative error between X and the solution is at ! most TOL. TOL should be nonnegative. ! ! Output, integer ( kind = 4 ) INFO, error flag. 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: ! 0, improper input parameters. ! 1, algorithm estimates that the relative error in the sum of squares ! is at most TOL. ! 2, algorithm estimates that the relative error between X and the ! solution is at most TOL. ! 3, conditions for INFO = 1 and INFO = 2 both hold. ! 4, FVEC is orthogonal to the columns of the jacobian to machine precision. ! 5, number of calls to FCN has reached or exceeded 200*(N+1). ! 6, TOL is too small. No further reduction in the sum of squares ! is possible. ! 7, TOL is too small. No further improvement in the approximate ! solution X is possible. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) diag(n) real ( kind = 8 ) epsfcn real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(m,n) real ( kind = 8 ) ftol real ( kind = 8 ) fvec(m) real ( kind = 8 ) gtol integer ( kind = 4 ) info integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) ldfjac integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) nprint real ( kind = 8 ) qtf(n) real ( kind = 8 ) tol real ( kind = 8 ) x(n) real ( kind = 8 ) xtol info = 0 if ( n <= 0 ) then return else if ( m < n ) then return else if ( tol < 0.0D+00 ) then return end if factor = 100.0D+00 maxfev = 200 * ( n + 1 ) ftol = tol xtol = tol gtol = 0.0D+00 epsfcn = 0.0D+00 mode = 1 nprint = 0 ldfjac = m call lmdif ( fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, & diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf ) if ( info == 8 ) then info = 4 end if return endsubroutine lmdif1 subroutine lmpar ( n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag ) !*****************************************************************************80 ! !! LMPAR computes a parameter for the Levenberg-Marquardt method. ! ! Discussion: ! ! Given an M by N matrix A, an N by N nonsingular diagonal ! matrix D, an M-vector B, and a positive number DELTA, ! the problem is to determine a value for the parameter ! PAR such that if X solves the system ! ! A*X = B, ! sqrt ( PAR ) * D * X = 0, ! ! in the least squares sense, and DXNORM is the euclidean ! norm of D*X, then either PAR is zero and ! ! ( DXNORM - DELTA ) <= 0.1 * DELTA, ! ! or PAR is positive and ! ! abs ( DXNORM - DELTA) <= 0.1 * DELTA. ! ! This function completes the solution of the problem ! if it is provided with the necessary information from the ! QR factorization, with column pivoting, of A. That is, if ! A*P = Q*R, where P is a permutation matrix, Q has orthogonal ! columns, and R is an upper triangular matrix with diagonal ! elements of nonincreasing magnitude, then LMPAR expects ! the full upper triangle of R, the permutation matrix P, ! and the first N components of Q'*B. On output ! LMPAR also provides an upper triangular matrix S such that ! ! P' * ( A' * A + PAR * D * D ) * P = S'* S. ! ! S is employed within LMPAR and may be of separate interest. ! ! Only a few iterations are generally needed for convergence ! of the algorithm. ! ! If, however, the limit of 10 iterations is reached, then the output ! PAR will contain the best value obtained so far. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 24 January 2014 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of R. ! ! Input/output, real ( kind = 8 ) R(LDR,N),the N by N matrix. The full ! upper triangle must contain the full upper triangle of the matrix R. ! On output the full upper triangle is unaltered, and the strict lower ! triangle contains the strict upper triangle (transposed) of the upper ! triangular matrix S. ! ! Input, integer ( kind = 4 ) LDR, the leading dimension of R. LDR must be ! no less than N. ! ! Input, integer ( kind = 4 ) IPVT(N), defines the permutation matrix P ! such that A*P = Q*R. Column J of P is column IPVT(J) of the ! identity matrix. ! ! Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D. ! ! Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'*B. ! ! Input, real ( kind = 8 ) DELTA, an upper bound on the euclidean norm ! of D*X. DELTA should be positive. ! ! Input/output, real ( kind = 8 ) PAR. On input an initial estimate of the ! Levenberg-Marquardt parameter. On output the final estimate. ! PAR should be nonnegative. ! ! Output, real ( kind = 8 ) X(N), the least squares solution of the system ! A*X = B, sqrt(PAR)*D*X = 0, for the output value of PAR. ! ! Output, real ( kind = 8 ) SDIAG(N), the diagonal elements of the upper ! triangular matrix S. ! implicit none integer ( kind = 4 ) ldr integer ( kind = 4 ) n real ( kind = 8 ) delta real ( kind = 8 ) diag(n) real ( kind = 8 ) dwarf real ( kind = 8 ) dxnorm !~ real ( kind = 8 ) enorm real ( kind = 8 ) gnorm real ( kind = 8 ) fp integer ( kind = 4 ) i integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) iter integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) nsing real ( kind = 8 ) par real ( kind = 8 ) parc real ( kind = 8 ) parl real ( kind = 8 ) paru real ( kind = 8 ) qnorm real ( kind = 8 ) qtb(n) real ( kind = 8 ) r(ldr,n) real ( kind = 8 ) sdiag(n) real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) x(n) ! ! DWARF is the smallest positive magnitude. ! dwarf = tiny ( dwarf ) ! ! Compute and store in X the Gauss-Newton direction. ! ! If the jacobian is rank-deficient, obtain a least squares solution. ! nsing = n do j = 1, n wa1(j) = qtb(j) if ( r(j,j) == 0.0D+00 .and. nsing == n ) then nsing = j - 1 end if if ( nsing < n ) then wa1(j) = 0.0D+00 end if end do do k = 1, nsing j = nsing - k + 1 wa1(j) = wa1(j) / r(j,j) temp = wa1(j) wa1(1:j-1) = wa1(1:j-1) - r(1:j-1,j) * temp end do do j = 1, n l = ipvt(j) x(l) = wa1(j) end do ! ! Initialize the iteration counter. ! Evaluate the function at the origin, and test ! for acceptance of the Gauss-Newton direction. ! iter = 0 wa2(1:n) = diag(1:n) * x(1:n) dxnorm = enorm ( n, wa2 ) fp = dxnorm - delta if ( fp <= 0.1D+00 * delta ) then if ( iter == 0 ) then par = 0.0D+00 end if return end if ! ! If the jacobian is not rank deficient, the Newton ! step provides a lower bound, PARL, for the zero of ! the function. ! ! Otherwise set this bound to zero. ! parl = 0.0D+00 if ( n <= nsing ) then do j = 1, n l = ipvt(j) wa1(j) = diag(l) * ( wa2(l) / dxnorm ) end do do j = 1, n sum2 = dot_product ( wa1(1:j-1), r(1:j-1,j) ) wa1(j) = ( wa1(j) - sum2 ) / r(j,j) end do temp = enorm ( n, wa1 ) parl = ( ( fp / delta ) / temp ) / temp end if ! ! Calculate an upper bound, PARU, for the zero of the function. ! do j = 1, n sum2 = dot_product ( qtb(1:j), r(1:j,j) ) l = ipvt(j) wa1(j) = sum2 / diag(l) end do gnorm = enorm ( n, wa1 ) paru = gnorm / delta if ( paru == 0.0D+00 ) then paru = dwarf / min ( delta, 0.1D+00 ) end if ! ! If the input PAR lies outside of the interval (PARL, PARU), ! set PAR to the closer endpoint. ! par = max ( par, parl ) par = min ( par, paru ) if ( par == 0.0D+00 ) then par = gnorm / dxnorm end if ! ! Beginning of an iteration. ! do iter = iter + 1 ! ! Evaluate the function at the current value of PAR. ! if ( par == 0.0D+00 ) then par = max ( dwarf, 0.001D+00 * paru ) end if wa1(1:n) = sqrt ( par ) * diag(1:n) call qrsolv ( n, r, ldr, ipvt, wa1, qtb, x, sdiag ) wa2(1:n) = diag(1:n) * x(1:n) dxnorm = enorm ( n, wa2 ) temp = fp fp = dxnorm - delta ! ! If the function is small enough, accept the current value of PAR. ! if ( abs ( fp ) <= 0.1D+00 * delta ) then exit end if ! ! Test for the exceptional cases where PARL ! is zero or the number of iterations has reached 10. ! if ( parl == 0.0D+00 .and. fp <= temp .and. temp < 0.0D+00 ) then exit else if ( iter == 10 ) then exit end if ! ! Compute the Newton correction. ! do j = 1, n l = ipvt(j) wa1(j) = diag(l) * ( wa2(l) / dxnorm ) end do do j = 1, n wa1(j) = wa1(j) / sdiag(j) temp = wa1(j) wa1(j+1:n) = wa1(j+1:n) - r(j+1:n,j) * temp end do temp = enorm ( n, wa1 ) parc = ( ( fp / delta ) / temp ) / temp ! ! Depending on the sign of the function, update PARL or PARU. ! if ( 0.0D+00 < fp ) then parl = max ( parl, par ) else if ( fp < 0.0D+00 ) then paru = min ( paru, par ) end if ! ! Compute an improved estimate for PAR. ! par = max ( parl, par + parc ) ! ! End of an iteration. ! end do ! ! Termination. ! if ( iter == 0 ) then par = 0.0D+00 end if return endsubroutine lmpar subroutine lmstr ( fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf ) !*****************************************************************************80 ! !! LMSTR minimizes M functions in N variables using Levenberg-Marquardt method. ! ! Discussion: ! ! LMSTR minimizes the sum of the squares of M nonlinear functions in ! N variables by a modification of the Levenberg-Marquardt algorithm ! which uses minimal storage. ! ! The user must provide a subroutine which calculates the functions and ! the rows of the jacobian. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions and the rows of the jacobian. ! FCN should have the form: ! subroutine fcn ( m, n, x, fvec, fjrow, iflag ) ! integer ( kind = 4 ) m ! integer ( kind = 4 ) n ! real ( kind = 8 ) fjrow(n) ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! If the input value of IFLAG is I > 1, calculate the (I-1)-st row of ! the jacobian at X, and return this vector in FJROW. ! To terminate the algorithm, set the output value of IFLAG negative. ! ! Input, integer ( kind = 4 ) M, the number of functions. ! ! Input, integer ( kind = 4 ) N, the number of variables. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array. The upper ! triangle of FJAC contains an upper triangular matrix R such that ! ! P' * ( JAC' * JAC ) * P = R' * R, ! ! where P is a permutation matrix and JAC is the final calculated jacobian. ! Column J of P is column IPVT(J) of the identity matrix. The lower ! triangular part of FJAC contains information generated during ! the computation of R. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least N. ! ! Input, real ( kind = 8 ) FTOL. Termination occurs when both the actual and ! predicted relative reductions in the sum of squares are at most FTOL. ! Therefore, FTOL measures the relative error desired in the sum of ! squares. FTOL should be nonnegative. ! ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error ! between two consecutive iterates is at most XTOL. XTOL should be ! nonnegative. ! ! Input, real ( kind = 8 ) GTOL. termination occurs when the cosine of the ! angle between FVEC and any column of the jacobian is at most GTOL in ! absolute value. Therefore, GTOL measures the orthogonality desired ! between the function vector and the columns of the jacobian. GTOL should ! be nonnegative. ! ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number ! of calls to FCN with IFLAG = 1 is at least MAXFEV by the end of ! an iteration. ! ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set ! internally. If MODE = 2, then DIAG must contain positive entries that ! serve as multiplicative scale factors for the variables. ! ! Input, integer ( kind = 4 ) MODE, scaling option. ! 1, variables will be scaled internally. ! 2, scaling is specified by the input DIAG vector. ! ! Input, real ( kind = 8 ) FACTOR, determines 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 (0.1, 100) with 100 the recommended value. ! ! Input, integer ( kind = 4 ) NPRINT, 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. ! ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated ! execution, INFO is set to the (negative) value of IFLAG. See the ! description of FCN. Otherwise, INFO is set as follows: ! 0, improper input parameters. ! 1, both actual and predicted relative reductions in the sum of squares ! are at most FTOL. ! 2, relative error between two consecutive iterates is at most XTOL. ! 3, conditions for INFO = 1 and INFO = 2 both hold. ! 4, the cosine of the angle between FVEC and any column of the jacobian ! is at most GTOL in absolute value. ! 5, number of calls to FCN with IFLAG = 1 has reached MAXFEV. ! 6, FTOL is too small. No further reduction in the sum of squares is ! possible. ! 7, XTOL is too small. No further improvement in the approximate ! solution X is possible. ! 8, GTOL is too small. FVEC is orthogonal to the columns of the ! jacobian to machine precision. ! ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN ! with IFLAG = 1. ! ! Output, integer ( kind = 4 ) NJEV, the number of calls to FCN ! with IFLAG = 2. ! ! Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P such ! that JAC * P = Q * R, where JAC is the final calculated jacobian, Q is ! orthogonal (not stored), and R is upper triangular. ! Column J of P is column IPVT(J) of the identity matrix. ! ! Output, real ( kind = 8 ) QTF(N), contains the first N elements of Q'*FVEC. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) actred real ( kind = 8 ) delta real ( kind = 8 ) diag(n) real ( kind = 8 ) dirder !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsmch real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) fnorm real ( kind = 8 ) fnorm1 real ( kind = 8 ) ftol real ( kind = 8 ) fvec(m) real ( kind = 8 ) gnorm real ( kind = 8 ) gtol integer ( kind = 4 ) i integer ( kind = 4 ) iflag integer ( kind = 4 ) info integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) iter integer ( kind = 4 ) j integer ( kind = 4 ) l integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) njev integer ( kind = 4 ) nprint real ( kind = 8 ) par logical pivot real ( kind = 8 ) pnorm real ( kind = 8 ) prered real ( kind = 8 ) qtf(n) real ( kind = 8 ) ratio logical sing real ( kind = 8 ) sum2 real ( kind = 8 ) temp real ( kind = 8 ) temp1 real ( kind = 8 ) temp2 real ( kind = 8 ) wa1(n) real ( kind = 8 ) wa2(n) real ( kind = 8 ) wa3(n) real ( kind = 8 ) wa4(m) real ( kind = 8 ) x(n) real ( kind = 8 ) xnorm real ( kind = 8 ) xtol epsmch = epsilon ( epsmch ) info = 0 iflag = 0 nfev = 0 njev = 0 ! ! Check the input parameters for errors. ! if ( n <= 0 ) then go to 340 else if ( m < n ) then go to 340 else if ( ldfjac < n ) then go to 340 else if ( ftol < 0.0D+00 ) then go to 340 else if ( xtol < 0.0D+00 ) then go to 340 else if ( gtol < 0.0D+00 ) then go to 340 else if ( maxfev <= 0 ) then go to 340 else if ( factor <= 0.0D+00 ) then go to 340 end if if ( mode == 2 ) then do j = 1, n if ( diag(j) <= 0.0D+00 ) then go to 340 end if end do end if ! ! Evaluate the function at the starting point and calculate its norm. ! iflag = 1 call fcn ( m, n, x, fvec, wa3, iflag ) nfev = 1 if ( iflag < 0 ) then go to 340 end if fnorm = enorm ( m, fvec ) ! ! Initialize Levenberg-Marquardt parameter and iteration counter. ! par = 0.0D+00 iter = 1 ! ! Beginning of the outer loop. ! 30 continue ! ! If requested, call FCN to enable printing of iterates. ! if ( 0 < nprint ) then iflag = 0 if ( mod ( iter-1, nprint ) == 0 ) then call fcn ( m, n, x, fvec, wa3, iflag ) end if if ( iflag < 0 ) then go to 340 end if end if ! ! Compute the QR factorization of the jacobian matrix calculated one row ! at a time, while simultaneously forming Q'* FVEC and storing ! the first N components in QTF. ! qtf(1:n) = 0.0D+00 fjac(1:n,1:n) = 0.0D+00 iflag = 2 do i = 1, m call fcn ( m, n, x, fvec, wa3, iflag ) if ( iflag < 0 ) then go to 340 end if temp = fvec(i) call rwupdt ( n, fjac, ldfjac, wa3, qtf, temp, wa1, wa2 ) iflag = iflag + 1 end do njev = njev + 1 ! ! If the jacobian is rank deficient, call QRFAC to ! reorder its columns and update the components of QTF. ! sing = .false. do j = 1, n if ( fjac(j,j) == 0.0D+00 ) then sing = .true. end if ipvt(j) = j wa2(j) = enorm ( j, fjac(1,j) ) end do if ( sing ) then pivot = .true. call qrfac ( n, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2 ) do j = 1, n if ( fjac(j,j) /= 0.0D+00 ) then sum2 = dot_product ( qtf(j:n), fjac(j:n,j) ) temp = - sum2 / fjac(j,j) qtf(j:n) = qtf(j:n) + fjac(j:n,j) * temp end if fjac(j,j) = wa1(j) end do end if ! ! On the first iteration ! if mode is 1, ! scale according to the norms of the columns of the initial jacobian. ! calculate the norm of the scaled X, ! initialize the step bound delta. ! if ( iter == 1 ) then if ( mode /= 2 ) then diag(1:n) = wa2(1:n) do j = 1, n if ( wa2(j) == 0.0D+00 ) then diag(j) = 1.0D+00 end if end do end if wa3(1:n) = diag(1:n) * x(1:n) xnorm = enorm ( n, wa3 ) delta = factor * xnorm if ( delta == 0.0D+00 ) then delta = factor end if end if ! ! Compute the norm of the scaled gradient. ! gnorm = 0.0D+00 if ( fnorm /= 0.0D+00 ) then do j = 1, n l = ipvt(j) if ( wa2(l) /= 0.0D+00 ) then sum2 = dot_product ( qtf(1:j), fjac(1:j,j) ) / fnorm gnorm = max ( gnorm, abs ( sum2 / wa2(l) ) ) end if end do end if ! ! Test for convergence of the gradient norm. ! if ( gnorm <= gtol ) then info = 4 go to 340 end if ! ! Rescale if necessary. ! if ( mode /= 2 ) then do j = 1, n diag(j) = max ( diag(j), wa2(j) ) end do end if ! ! Beginning of the inner loop. ! 240 continue ! ! Determine the Levenberg-Marquardt parameter. ! call lmpar ( n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2 ) ! ! Store the direction P and X + P. ! Calculate the norm of P. ! wa1(1:n) = -wa1(1:n) wa2(1:n) = x(1:n) + wa1(1:n) wa3(1:n) = diag(1:n) * wa1(1:n) pnorm = enorm ( n, wa3 ) ! ! On the first iteration, adjust the initial step bound. ! if ( iter == 1 ) then delta = min ( delta, pnorm ) end if ! ! Evaluate the function at X + P and calculate its norm. ! iflag = 1 call fcn ( m, n, wa2, wa4, wa3, iflag ) nfev = nfev + 1 if ( iflag < 0 ) then go to 340 end if fnorm1 = enorm ( m, wa4 ) ! ! Compute the scaled actual reduction. ! if ( 0.1D+00 * fnorm1 < fnorm ) then actred = 1.0D+00 - ( fnorm1 / fnorm ) ** 2 else actred = -1.0D+00 end if ! ! Compute the scaled predicted reduction and ! the scaled directional derivative. ! do j = 1, n wa3(j) = 0.0D+00 l = ipvt(j) temp = wa1(l) wa3(1:j) = wa3(1:j) + fjac(1:j,j) * temp end do temp1 = enorm ( n, wa3 ) / fnorm temp2 = ( sqrt(par) * pnorm ) / fnorm prered = temp1 ** 2 + temp2 ** 2 / 0.5D+00 dirder = - ( temp1 ** 2 + temp2 ** 2 ) ! ! Compute the ratio of the actual to the predicted reduction. ! if ( prered /= 0.0D+00 ) then ratio = actred / prered else ratio = 0.0D+00 end if ! ! Update the step bound. ! if ( ratio <= 0.25D+00 ) then if ( actred >= 0.0D+00 ) then temp = 0.5D+00 else temp = 0.5D+00 * dirder / ( dirder + 0.5D+00 * actred ) end if if ( 0.1D+00 * fnorm1 >= fnorm .or. temp < 0.1D+00 ) then temp = 0.1D+00 end if delta = temp * min ( delta, pnorm / 0.1D+00 ) par = par / temp else if ( par == 0.0D+00 .or. ratio >= 0.75D+00 ) then delta = pnorm / 0.5D+00 par = 0.5D+00 * par end if end if ! ! Test for successful iteration. ! if ( ratio >= 0.0001D+00 ) then x(1:n) = wa2(1:n) wa2(1:n) = diag(1:n) * x(1:n) fvec(1:m) = wa4(1:m) xnorm = enorm ( n, wa2 ) fnorm = fnorm1 iter = iter + 1 end if ! ! Tests for convergence, termination and stringent tolerances. ! if ( abs ( actred ) <= ftol .and. prered <= ftol & .and. 0.5D+00 * ratio <= 1.0D+00 ) then info = 1 end if if ( delta <= xtol * xnorm ) then info = 2 end if if ( abs ( actred ) <= ftol .and. prered <= ftol & .and. 0.5D+00 * ratio <= 1.0D+00 .and. info == 2 ) then info = 3 end if if ( info /= 0 ) then go to 340 end if if ( nfev >= maxfev) then info = 5 end if if ( abs ( actred ) <= epsmch .and. prered <= epsmch & .and. 0.5D+00 * ratio <= 1.0D+00 ) then info = 6 end if if ( delta <= epsmch * xnorm ) then info = 7 end if if ( gnorm <= epsmch ) then info = 8 end if if ( info /= 0 ) then go to 340 end if ! ! End of the inner loop. Repeat if iteration unsuccessful. ! if ( ratio < 0.0001D+00 ) then go to 240 end if ! ! End of the outer loop. ! go to 30 340 continue ! ! Termination, either normal or user imposed. ! if ( iflag < 0 ) then info = iflag end if iflag = 0 if ( 0 < nprint ) then call fcn ( m, n, x, fvec, wa3, iflag ) end if return endsubroutine lmstr subroutine lmstr1 ( fcn, m, n, x, fvec, fjac, ldfjac, tol, info ) !*****************************************************************************80 ! !! LMSTR1 minimizes M functions in N variables using Levenberg-Marquardt method. ! ! Discussion: ! ! LMSTR1 minimizes the sum of the squares of M nonlinear functions in ! N variables by a modification of the Levenberg-Marquardt algorithm ! which uses minimal storage. ! ! This is done by using the more general least-squares solver ! LMSTR. The user must provide a subroutine which calculates ! the functions and the rows of the jacobian. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 19 August 2016 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, external FCN, the name of the user-supplied subroutine which ! calculates the functions and the rows of the jacobian. ! FCN should have the form: ! subroutine fcn ( m, n, x, fvec, fjrow, iflag ) ! integer ( kind = 4 ) m ! integer ( kind = 4 ) n ! real ( kind = 8 ) fjrow(n) ! real ( kind = 8 ) fvec(m) ! integer ( kind = 4 ) iflag ! real ( kind = 8 ) x(n) ! If IFLAG = 0 on input, then FCN is only being called to allow the user ! to print out the current iterate. ! If IFLAG = 1 on input, FCN should calculate the functions at X and ! return this vector in FVEC. ! If the input value of IFLAG is I > 1, calculate the (I-1)-st row of ! the jacobian at X, and return this vector in FJROW. ! To terminate the algorithm, set the output value of IFLAG negative. ! ! Input, integer ( kind = 4 ) M, the number of functions. ! ! Input, integer ( kind = 4 ) N, the number of variables. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial ! estimate of the solution vector. On output X contains the final ! estimate of the solution vector. ! ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X. ! ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array. The upper ! triangle contains an upper triangular matrix R such that ! ! P' * ( JAC' * JAC ) * P = R' * R, ! ! where P is a permutation matrix and JAC is the final calculated ! jacobian. Column J of P is column IPVT(J) of the identity matrix. ! The lower triangular part of FJAC contains information generated ! during the computation of R. ! ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC. ! LDFJAC must be at least N. ! ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm ! estimates either that the relative error in the sum of squares is at ! most TOL or that the relative error between X and the solution is at ! most TOL. TOL should be nonnegative. ! ! Output, integer ( kind = 4 ) INFO, error flag. 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: ! 0, improper input parameters. ! 1, algorithm estimates that the relative error in the sum of squares ! is at most TOL. ! 2, algorithm estimates that the relative error between X and the ! solution is at most TOL. ! 3, conditions for INFO = 1 and INFO = 2 both hold. ! 4, FVEC is orthogonal to the columns of the jacobian to machine precision. ! 5, number of calls to FCN with IFLAG = 1 has reached 100*(N+1). ! 6, TOL is too small. No further reduction in the sum of squares ! is possible. ! 7, TOL is too small. No further improvement in the approximate ! solution X is possible. ! implicit none integer ( kind = 4 ) ldfjac integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) diag(n) real ( kind = 8 ) factor external fcn real ( kind = 8 ) fjac(ldfjac,n) real ( kind = 8 ) ftol real ( kind = 8 ) fvec(m) real ( kind = 8 ) gtol integer ( kind = 4 ) info integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) maxfev integer ( kind = 4 ) mode integer ( kind = 4 ) nfev integer ( kind = 4 ) njev integer ( kind = 4 ) nprint real ( kind = 8 ) qtf(n) real ( kind = 8 ) tol real ( kind = 8 ) x(n) real ( kind = 8 ) xtol if ( n <= 0 ) then info = 0 return end if if ( m < n ) then info = 0 return end if if ( ldfjac < n ) then info = 0 return end if if ( tol < 0.0D+00 ) then info = 0 return end if fvec(1:n) = 0.0D+00 fjac(1:ldfjac,1:n) = 0.0D+00 ftol = tol xtol = tol gtol = 0.0D+00 maxfev = 100 * ( n + 1 ) diag(1:n) = 0.0D+00 mode = 1 factor = 100.0D+00 nprint = 0 info = 0 nfev = 0 njev = 0 ipvt(1:n) = 0 qtf(1:n) = 0.0D+00 call lmstr ( fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf ) if ( info == 8 ) then info = 4 end if return endsubroutine lmstr1 subroutine qform ( m, n, q, ldq ) !*****************************************************************************80 ! !! QFORM produces the explicit QR factorization of a matrix. ! ! Discussion: ! ! The QR factorization of a matrix is usually accumulated in implicit ! form, that is, as a series of orthogonal transformations of the ! original matrix. This routine carries out those transformations, ! to explicitly exhibit the factorization constructed by QRFAC. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, is a positive integer input variable set ! to the number of rows of A and the order of Q. ! ! Input, integer ( kind = 4 ) N, is a positive integer input variable set ! to the number of columns of A. ! ! Input/output, real ( kind = 8 ) Q(LDQ,M). Q is an M by M array. ! On input the full lower trapezoid in the first min(M,N) columns of Q ! contains the factored form. ! On output, Q has been accumulated into a square matrix. ! ! Input, integer ( kind = 4 ) LDQ, is a positive integer input variable ! not less than M which specifies the leading dimension of the array Q. ! implicit none integer ( kind = 4 ) ldq integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) minmn real ( kind = 8 ) q(ldq,m) real ( kind = 8 ) temp real ( kind = 8 ) wa(m) minmn = min ( m, n ) do j = 2, minmn q(1:j-1,j) = 0.0D+00 end do ! ! Initialize remaining columns to those of the identity matrix. ! q(1:m,n+1:m) = 0.0D+00 do j = n + 1, m q(j,j) = 1.0D+00 end do ! ! Accumulate Q from its factored form. ! do l = 1, minmn k = minmn - l + 1 wa(k:m) = q(k:m,k) q(k:m,k) = 0.0D+00 q(k,k) = 1.0D+00 if ( wa(k) /= 0.0D+00 ) then do j = k, m temp = dot_product ( wa(k:m), q(k:m,j) ) / wa(k) q(k:m,j) = q(k:m,j) - temp * wa(k:m) end do end if end do return endsubroutine qform subroutine qrfac ( m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm ) !*****************************************************************************80 ! !! QRFAC computes a QR factorization using Householder transformations. ! ! Discussion: ! ! This function uses Householder transformations with optional column ! pivoting to compute a QR factorization of the ! M by N matrix A. That is, QRFAC determines an orthogonal ! matrix Q, a permutation matrix P, and an upper trapezoidal ! matrix R with diagonal elements of nonincreasing magnitude, ! such that A*P = Q*R. ! ! The Householder transformation for column K, K = 1,2,...,min(M,N), ! is of the form ! ! I - ( 1 / U(K) ) * U * U' ! ! where U has zeros in the first K-1 positions. ! ! The form of this transformation and the method of pivoting first ! appeared in the corresponding LINPACK routine. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows of A. ! ! Input, integer ( kind = 4 ) N, the number of columns of A. ! ! Input/output, real ( kind = 8 ) A(LDA,N), the M by N array. ! On input, A contains the matrix for which the QR factorization is to ! be computed. On output, the strict upper trapezoidal part of A contains ! the strict upper trapezoidal part of R, and the lower trapezoidal ! part of A contains a factored form of Q, the non-trivial elements of ! the U vectors described above. ! ! Input, integer ( kind = 4 ) LDA, the leading dimension of A, which must ! be no less than M. ! ! Input, logical PIVOT, is TRUE if column pivoting is to be carried out. ! ! Output, integer ( kind = 4 ) IPVT(LIPVT), defines the permutation matrix P ! such that A*P = Q*R. Column J of P is column IPVT(J) of the identity ! matrix. If PIVOT is false, IPVT is not referenced. ! ! Input, integer ( kind = 4 ) LIPVT, the dimension of IPVT, which should ! be N if pivoting is used. ! ! Output, real ( kind = 8 ) RDIAG(N), contains the diagonal elements of R. ! ! Output, real ( kind = 8 ) ACNORM(N), the norms of the corresponding ! columns of the input matrix A. If this information is not needed, ! then ACNORM can coincide with RDIAG. ! implicit none integer ( kind = 4 ) lda integer ( kind = 4 ) lipvt integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(lda,n) real ( kind = 8 ) acnorm(n) real ( kind = 8 ) ajnorm !~ real ( kind = 8 ) enorm real ( kind = 8 ) epsmch integer ( kind = 4 ) i integer ( kind = 4 ) i4_temp integer ( kind = 4 ) ipvt(lipvt) integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) kmax integer ( kind = 4 ) minmn logical pivot real ( kind = 8 ) r8_temp(m) real ( kind = 8 ) rdiag(n) real ( kind = 8 ) temp real ( kind = 8 ) wa(n) epsmch = epsilon ( epsmch ) ! ! Compute the initial column norms and initialize several arrays. ! do j = 1, n acnorm(j) = enorm( m, a(1:m,j) ) end do rdiag(1:n) = acnorm(1:n) wa(1:n) = acnorm(1:n) if ( pivot ) then do j = 1, n ipvt(j) = j end do end if ! ! Reduce A to R with Householder transformations. ! minmn = min ( m, n ) do j = 1, minmn ! ! Bring the column of largest norm into the pivot position. ! if ( pivot ) then kmax = j do k = j, n if ( rdiag(kmax) < rdiag(k) ) then kmax = k end if end do if ( kmax /= j ) then r8_temp(1:m) = a(1:m,j) a(1:m,j) = a(1:m,kmax) a(1:m,kmax) = r8_temp(1:m) rdiag(kmax) = rdiag(j) wa(kmax) = wa(j) i4_temp = ipvt(j) ipvt(j) = ipvt(kmax) ipvt(kmax) = i4_temp end if end if ! ! Compute the Householder transformation to reduce the ! J-th column of A to a multiple of the J-th unit vector. ! ajnorm = enorm ( m-j+1, a(j,j) ) if ( ajnorm /= 0.0D+00 ) then if ( a(j,j) < 0.0D+00 ) then ajnorm = -ajnorm end if a(j:m,j) = a(j:m,j) / ajnorm a(j,j) = a(j,j) + 1.0D+00 ! ! Apply the transformation to the remaining columns and update the norms. ! do k = j + 1, n temp = dot_product ( a(j:m,j), a(j:m,k) ) / a(j,j) a(j:m,k) = a(j:m,k) - temp * a(j:m,j) if ( pivot .and. rdiag(k) /= 0.0D+00 ) then temp = a(j,k) / rdiag(k) rdiag(k) = rdiag(k) * sqrt ( max ( 0.0D+00, 1.0D+00-temp ** 2 ) ) if ( 0.05D+00 * ( rdiag(k) / wa(k) ) ** 2 <= epsmch ) then rdiag(k) = enorm ( m-j, a(j+1,k) ) wa(k) = rdiag(k) end if end if end do end if rdiag(j) = - ajnorm end do return endsubroutine qrfac subroutine qrsolv ( n, r, ldr, ipvt, diag, qtb, x, sdiag ) !*****************************************************************************80 ! !! QRSOLV solves a rectangular linear system A*x=b in the least squares sense. ! ! Discussion: ! ! Given an M by N matrix A, an N by N diagonal matrix D, ! and an M-vector B, the problem is to determine an X which ! solves the system ! ! A*X = B ! D*X = 0 ! ! in the least squares sense. ! ! This function completes the solution of the problem ! if it is provided with the necessary information from the ! QR factorization, with column pivoting, of A. That is, if ! Q*P = Q*R, where P is a permutation matrix, Q has orthogonal ! columns, and R is an upper triangular matrix with diagonal ! elements of nonincreasing magnitude, then QRSOLV expects ! the full upper triangle of R, the permutation matrix p, ! and the first N components of Q'*B. ! ! The system is then equivalent to ! ! R*Z = Q'*B ! P'*D*P*Z = 0 ! ! where X = P*Z. If this system does not have full rank, ! then a least squares solution is obtained. On output QRSOLV ! also provides an upper triangular matrix S such that ! ! P'*(A'*A + D*D)*P = S'*S. ! ! S is computed within QRSOLV and may be of separate interest. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of R. ! ! Input/output, real ( kind = 8 ) R(LDR,N), the N by N matrix. ! On input the full upper triangle must contain the full upper triangle ! of the matrix R. On output the full upper triangle is unaltered, and ! the strict lower triangle contains the strict upper triangle ! (transposed) of the upper triangular matrix S. ! ! Input, integer ( kind = 4 ) LDR, the leading dimension of R, which must be ! at least N. ! ! Input, integer ( kind = 4 ) IPVT(N), defines the permutation matrix P such ! that A*P = Q*R. Column J of P is column IPVT(J) of the identity matrix. ! ! Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D. ! ! Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'*B. ! ! Output, real ( kind = 8 ) X(N), the least squares solution. ! ! Output, real ( kind = 8 ) SDIAG(N), the diagonal elements of the upper ! triangular matrix S. ! implicit none integer ( kind = 4 ) ldr integer ( kind = 4 ) n real ( kind = 8 ) c real ( kind = 8 ) cotan real ( kind = 8 ) diag(n) integer ( kind = 4 ) i integer ( kind = 4 ) ipvt(n) integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) nsing real ( kind = 8 ) qtb(n) real ( kind = 8 ) qtbpj real ( kind = 8 ) r(ldr,n) real ( kind = 8 ) s real ( kind = 8 ) sdiag(n) real ( kind = 8 ) sum2 real ( kind = 8 ) t real ( kind = 8 ) temp real ( kind = 8 ) wa(n) real ( kind = 8 ) x(n) ! ! Copy R and Q'*B to preserve input and initialize S. ! ! In particular, save the diagonal elements of R in X. ! do j = 1, n r(j:n,j) = r(j,j:n) x(j) = r(j,j) end do wa(1:n) = qtb(1:n) ! ! Eliminate the diagonal matrix D using a Givens rotation. ! do j = 1, n ! ! Prepare the row of D to be eliminated, locating the ! diagonal element using P from the QR factorization. ! l = ipvt(j) if ( diag(l) /= 0.0D+00 ) then sdiag(j:n) = 0.0D+00 sdiag(j) = diag(l) ! ! The transformations to eliminate the row of D ! modify only a single element of Q'*B ! beyond the first N, which is initially zero. ! qtbpj = 0.0D+00 do k = j, n ! ! Determine a Givens rotation which eliminates the ! appropriate element in the current row of D. ! if ( sdiag(k) /= 0.0D+00 ) then if ( abs ( r(k,k) ) < abs ( sdiag(k) ) ) then cotan = r(k,k) / sdiag(k) s = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * cotan ** 2 ) c = s * cotan else t = sdiag(k) / r(k,k) c = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * t ** 2 ) s = c * t end if ! ! Compute the modified diagonal element of R and ! the modified element of (Q'*B,0). ! r(k,k) = c * r(k,k) + s * sdiag(k) temp = c * wa(k) + s * qtbpj qtbpj = - s * wa(k) + c * qtbpj wa(k) = temp ! ! Accumulate the tranformation in the row of S. ! do i = k+1, n temp = c * r(i,k) + s * sdiag(i) sdiag(i) = - s * r(i,k) + c * sdiag(i) r(i,k) = temp end do end if end do end if ! ! Store the diagonal element of S and restore ! the corresponding diagonal element of R. ! sdiag(j) = r(j,j) r(j,j) = x(j) end do ! ! Solve the triangular system for Z. If the system is ! singular, then obtain a least squares solution. ! nsing = n do j = 1, n if ( sdiag(j) == 0.0D+00 .and. nsing == n ) then nsing = j - 1 end if if ( nsing < n ) then wa(j) = 0.0D+00 end if end do do j = nsing, 1, -1 sum2 = dot_product ( wa(j+1:nsing), r(j+1:nsing,j) ) wa(j) = ( wa(j) - sum2 ) / sdiag(j) end do ! ! Permute the components of Z back to components of X. ! do j = 1, n l = ipvt(j) x(l) = wa(j) end do return endsubroutine qrsolv subroutine r1mpyq ( m, n, a, lda, v, w ) !*****************************************************************************80 ! !! R1MPYQ computes A*Q, where Q is the product of Householder transformations. ! ! Discussion: ! ! Given an M by N matrix A, this function computes A*Q where ! Q is the product of 2*(N - 1) transformations ! ! GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) ! ! and GV(I), GW(I) are Givens rotations in the (I,N) plane which ! eliminate elements in the I-th and N-th planes, respectively. ! Q itself is not given, rather the information to recover the ! GV, GW rotations is supplied. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows of A. ! ! Input, integer ( kind = 4 ) N, the number of columns of A. ! ! Input/output, real ( kind = 8 ) A(LDA,N), the M by N array. ! On input, the matrix A to be postmultiplied by the orthogonal matrix Q. ! On output, the value of A*Q. ! ! Input, integer ( kind = 4 ) LDA, the leading dimension of A, which must not ! be less than M. ! ! Input, real ( kind = 8 ) V(N), W(N), contain the information necessary ! to recover the Givens rotations GV and GW. ! implicit none integer ( kind = 4 ) lda integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(lda,n) real ( kind = 8 ) c integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) s real ( kind = 8 ) temp real ( kind = 8 ) v(n) real ( kind = 8 ) w(n) ! ! Apply the first set of Givens rotations to A. ! do j = n - 1, 1, -1 if ( 1.0D+00 < abs ( v(j) ) ) then c = 1.0D+00 / v(j) s = sqrt ( 1.0D+00 - c ** 2 ) else s = v(j) c = sqrt ( 1.0D+00 - s ** 2 ) end if do i = 1, m temp = c * a(i,j) - s * a(i,n) a(i,n) = s * a(i,j) + c * a(i,n) a(i,j) = temp end do end do ! ! Apply the second set of Givens rotations to A. ! do j = 1, n - 1 if ( abs ( w(j) ) > 1.0D+00 ) then c = 1.0D+00 / w(j) s = sqrt ( 1.0D+00 - c ** 2 ) else s = w(j) c = sqrt ( 1.0D+00 - s ** 2 ) end if do i = 1, m temp = c * a(i,j) + s * a(i,n) a(i,n) = - s * a(i,j) + c * a(i,n) a(i,j) = temp end do end do return endsubroutine r1mpyq subroutine r1updt ( m, n, s, ls, u, v, w, sing ) !*****************************************************************************80 ! !! R1UPDT re-triangularizes a matrix after a rank one update. ! ! Discussion: ! ! Given an M by N lower trapezoidal matrix S, an M-vector U, and an ! N-vector V, the problem is to determine an orthogonal matrix Q such that ! ! (S + U * V' ) * Q ! ! is again lower trapezoidal. ! ! This function determines Q as the product of 2 * (N - 1) ! transformations ! ! GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) ! ! where GV(I), GW(I) are Givens rotations in the (I,N) plane ! which eliminate elements in the I-th and N-th planes, ! respectively. Q itself is not accumulated, rather the ! information to recover the GV and GW rotations is returned. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows of S. ! ! Input, integer ( kind = 4 ) N, the number of columns of S. ! N must not exceed M. ! ! Input/output, real ( kind = 8 ) S(LS). On input, the lower trapezoidal ! matrix S stored by columns. On output S contains the lower trapezoidal ! matrix produced as described above. ! ! Input, integer ( kind = 4 ) LS, the length of the S array. LS must be at ! least (N*(2*M-N+1))/2. ! ! Input, real ( kind = 8 ) U(M), the U vector. ! ! Input/output, real ( kind = 8 ) V(N). On input, V must contain the ! vector V. On output V contains the information necessary to recover the ! Givens rotations GV described above. ! ! Output, real ( kind = 8 ) W(M), contains information necessary to ! recover the Givens rotations GW described above. ! ! Output, logical SING, is set to TRUE if any of the diagonal elements ! of the output S are zero. Otherwise SING is set FALSE. ! implicit none integer ( kind = 4 ) ls integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) cos real ( kind = 8 ) cotan real ( kind = 8 ) giant integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) l real ( kind = 8 ) s(ls) real ( kind = 8 ) sin logical sing real ( kind = 8 ) tan real ( kind = 8 ) tau real ( kind = 8 ) temp real ( kind = 8 ) u(m) real ( kind = 8 ) v(n) real ( kind = 8 ) w(m) ! ! GIANT is the largest magnitude. ! giant = huge ( giant ) ! ! Initialize the diagonal element pointer. ! jj = ( n * ( 2 * m - n + 1 ) ) / 2 - ( m - n ) ! ! Move the nontrivial part of the last column of S into W. ! l = jj do i = n, m w(i) = s(l) l = l + 1 end do ! ! Rotate the vector V into a multiple of the N-th unit vector ! in such a way that a spike is introduced into W. ! do j = n - 1, 1, -1 jj = jj - ( m - j + 1 ) w(j) = 0.0D+00 if ( v(j) /= 0.0D+00 ) then ! ! Determine a Givens rotation which eliminates the ! J-th element of V. ! if ( abs ( v(n) ) < abs ( v(j) ) ) then cotan = v(n) / v(j) sin = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * cotan ** 2 ) cos = sin * cotan tau = 1.0D+00 if ( abs ( cos ) * giant > 1.0D+00 ) then tau = 1.0D+00 / cos end if else tan = v(j) / v(n) cos = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * tan ** 2 ) sin = cos * tan tau = sin end if ! ! Apply the transformation to V and store the information ! necessary to recover the Givens rotation. ! v(n) = sin * v(j) + cos * v(n) v(j) = tau ! ! Apply the transformation to S and extend the spike in W. ! l = jj do i = j, m temp = cos * s(l) - sin * w(i) w(i) = sin * s(l) + cos * w(i) s(l) = temp l = l + 1 end do end if end do ! ! Add the spike from the rank 1 update to W. ! w(1:m) = w(1:m) + v(n) * u(1:m) ! ! Eliminate the spike. ! sing = .false. do j = 1, n-1 if ( w(j) /= 0.0D+00 ) then ! ! Determine a Givens rotation which eliminates the ! J-th element of the spike. ! if ( abs ( s(jj) ) < abs ( w(j) ) ) then cotan = s(jj) / w(j) sin = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * cotan ** 2 ) cos = sin * cotan if ( 1.0D+00 < abs ( cos ) * giant ) then tau = 1.0D+00 / cos else tau = 1.0D+00 end if else tan = w(j) / s(jj) cos = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * tan ** 2 ) sin = cos * tan tau = sin end if ! ! Apply the transformation to S and reduce the spike in W. ! l = jj do i = j, m temp = cos * s(l) + sin * w(i) w(i) = - sin * s(l) + cos * w(i) s(l) = temp l = l + 1 end do ! ! Store the information necessary to recover the Givens rotation. ! w(j) = tau end if ! ! Test for zero diagonal elements in the output S. ! if ( s(jj) == 0.0D+00 ) then sing = .true. end if jj = jj + ( m - j + 1 ) end do ! ! Move W back into the last column of the output S. ! l = jj do i = n, m s(l) = w(i) l = l + 1 end do if ( s(jj) == 0.0D+00 ) then sing = .true. end if return endsubroutine r1updt function r8_uniform_01 ( seed ) !*****************************************************************************80 ! !! R8_UNIFORM_01 returns a unit pseudorandom R8. ! ! Discussion: ! ! An R8 is a real ( kind = 8 ) value. ! ! For now, the input quantity SEED is an integer variable. ! ! This routine implements the recursion ! ! seed = 16807 * seed mod ( 2^31 - 1 ) ! r8_uniform_01 = seed / ( 2^31 - 1 ) ! ! The integer arithmetic never requires more than 32 bits, ! including a sign bit. ! ! If the initial seed is 12345, then the first three computations are ! ! Input Output R8_UNIFORM_01 ! SEED SEED ! ! 12345 207482415 0.096616 ! 207482415 1790989824 0.833995 ! 1790989824 2035175616 0.947702 ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 05 July 2006 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley Interscience, page 95, 1998. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R8_UNIFORM_01, a new pseudorandom variate, ! strictly between 0 and 1. ! implicit none integer ( kind = 4 ), parameter :: i4_huge = 2147483647 integer ( kind = 4 ) k real ( kind = 8 ) r8_uniform_01 integer ( kind = 4 ) seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r8_uniform_01 = real ( seed, kind = 8 ) * 4.656612875D-10 return endfunction r8_uniform_01 subroutine r8mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_PRINT prints an R8MAT. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 12 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows in A. ! ! Input, integer ( kind = 4 ) N, the number of columns in A. ! ! Input, real ( kind = 8 ) A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = * ) title call r8mat_print_some ( m, n, a, 1, 1, m, n, title ) return endsubroutine r8mat_print subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_PRINT_SOME prints some of an R8MAT. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return end if do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i8,6x)' ) j end do write ( *, '('' Col '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 if ( a(i,j) == real ( int ( a(i,j) ), kind = 8 ) ) then write ( ctemp(j2), '(f8.0,6x)' ) a(i,j) else write ( ctemp(j2), '(g14.6)' ) a(i,j) end if end do write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc ) end do end do return endsubroutine r8mat_print_some subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 22 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of components of the vector. ! ! Input, real ( kind = 8 ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,g16.8)' ) i, a(i) end do return endsubroutine r8vec_print subroutine rwupdt ( n, r, ldr, w, b, alpha, c, s ) !*****************************************************************************80 ! !! RWUPDT computes the decomposition of triangular matrix augmented by one row. ! ! Discussion: ! ! Given an N by N upper triangular matrix R, this function ! computes the QR decomposition of the matrix formed when a row ! is added to R. If the row is specified by the vector W, then ! RWUPDT determines an orthogonal matrix Q such that when the ! N+1 by N matrix composed of R augmented by W is premultiplied ! by Q', the resulting matrix is upper trapezoidal. ! The matrix Q' is the product of N transformations ! ! G(N)*G(N-1)* ... *G(1) ! ! where G(I) is a Givens rotation in the (I,N+1) plane which eliminates ! elements in the (N+1)-st plane. RWUPDT also computes the product ! Q'*C where C is the (N+1)-vector (B,ALPHA). ! ! Q itself is not accumulated, rather the information to recover the G ! rotations is supplied. ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 06 April 2010 ! ! Author: ! ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jorge More, Burton Garbow, Kenneth Hillstrom, ! User Guide for MINPACK-1, ! Technical Report ANL-80-74, ! Argonne National Laboratory, 1980. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of R. ! ! Input/output, real ( kind = 8 ) R(LDR,N). On input the upper triangular ! part of R must contain the matrix to be updated. On output R contains the ! updated triangular matrix. ! ! Input, integer ( kind = 4 ) LDR, the leading dimension of the array R. ! LDR must not be less than N. ! ! Input, real ( kind = 8 ) W(N), the row vector to be added to R. ! ! Input/output, real ( kind = 8 ) B(N). On input, the first N elements ! of the vector C. On output the first N elements of the vector Q'*C. ! ! Input/output, real ( kind = 8 ) ALPHA. On input, the (N+1)-st element ! of the vector C. On output the (N+1)-st element of the vector Q'*C. ! ! Output, real ( kind = 8 ) C(N), S(N), the cosines and sines of the ! transforming Givens rotations. ! implicit none integer ( kind = 4 ) ldr integer ( kind = 4 ) n real ( kind = 8 ) alpha real ( kind = 8 ) b(n) real ( kind = 8 ) c(n) real ( kind = 8 ) cotan integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) r(ldr,n) real ( kind = 8 ) rowj real ( kind = 8 ) s(n) real ( kind = 8 ) tan real ( kind = 8 ) temp real ( kind = 8 ) w(n) do j = 1, n rowj = w(j) ! ! Apply the previous transformations to R(I,J), I=1,2,...,J-1, and to W(J). ! do i = 1, j - 1 temp = c(i) * r(i,j) + s(i) * rowj rowj = - s(i) * r(i,j) + c(i) * rowj r(i,j) = temp end do ! ! Determine a Givens rotation which eliminates W(J). ! c(j) = 1.0D+00 s(j) = 0.0D+00 if ( rowj /= 0.0D+00 ) then if ( abs ( r(j,j) ) < abs ( rowj ) ) then cotan = r(j,j) / rowj s(j) = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * cotan ** 2 ) c(j) = s(j) * cotan else tan = rowj / r(j,j) c(j) = 0.5D+00 / sqrt ( 0.25D+00 + 0.25D+00 * tan ** 2 ) s(j) = c(j) * tan end if ! ! Apply the current transformation to R(J,J), B(J), and ALPHA. ! r(j,j) = c(j) * r(j,j) + s(j) * rowj temp = c(j) * b(j) + s(j) * alpha alpha = - s(j) * b(j) + c(j) * alpha b(j) = temp end if end do return endsubroutine rwupdt subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code may freely be copied, modified, and used for any purpose. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer ( kind = 4 ) d integer ( kind = 4 ) h integer ( kind = 4 ) m integer ( kind = 4 ) mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer ( kind = 4 ) n integer ( kind = 4 ) s integer ( kind = 4 ) values(8) integer ( kind = 4 ) y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return endsubroutine timestamp endmodule minpack