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 subroutine 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 ap = qr, 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 transpose)*b. on output lmpar also provides an upper triangular matrix s such that
t t t
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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n | a positive integer input variable set to the order of r. |
||
| real(kind=wp), | intent(inout) | :: | r(Ldr,n) | an n by n array. 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. |
||
| integer, | intent(in) | :: | Ldr | a positive integer input variable not less than n which specifies the leading dimension of the array r. |
||
| integer, | intent(in) | :: | Ipvt(n) | an integer input array of length n which defines the permutation matrix p such that ap = qr. column j of p is column ipvt(j) of the identity matrix. |
||
| real(kind=wp), | intent(in) | :: | Diag(n) | an input array of length n which must contain the diagonal elements of the matrix d. |
||
| real(kind=wp), | intent(in) | :: | Qtb(n) | an input array of length n which must contain the first n elements of the vector (q transpose)*b. |
||
| real(kind=wp) | :: | Delta | a positive input variable which specifies an upper bound on the euclidean norm of d*x. |
|||
| real(kind=wp), | intent(inout) | :: | Par | a nonnegative variable. on input par contains an initial estimate of the levenberg-marquardt parameter. on output par contains the final estimate. |
||
| real(kind=wp), | intent(out) | :: | x(n) | an output array of length n which contains the least squares solution of the system ax = b, sqrt(par)d*x = 0, for the output par. |
||
| real(kind=wp), | intent(out) | :: | Sdiag(n) | an output array of length n which contains the diagonal elements of the upper triangular matrix s. |
||
| real(kind=wp), | intent(inout) | :: | Wa1(n) | work array of length n. |
||
| real(kind=wp), | intent(inout) | :: | Wa2(n) | work array of length n. |
subroutine lmpar(n, r, Ldr, Ipvt, Diag, Qtb, Delta, Par, x, Sdiag, Wa1, Wa2)
implicit none
integer, intent(in) :: n !! a positive integer input variable set to the order of r.
integer, intent(in) :: Ldr !! a positive integer input variable not less than n
!! which specifies the leading dimension of the array r.
integer, intent(in) :: Ipvt(n) !! an integer input array of length n which defines the
!! permutation matrix p such that a*p = q*r. column j of p
!! is column ipvt(j) of the identity matrix.
real(wp) :: Delta !! a positive input variable which specifies an upper
!! bound on the euclidean norm of d*x.
real(wp), intent(inout) :: Par !! a nonnegative variable. on input par contains an
!! initial estimate of the levenberg-marquardt parameter.
!! on output par contains the final estimate.
real(wp), intent(inout) :: r(Ldr, n) !! an n by n array. 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.
real(wp), intent(in) :: Diag(n) !! an input array of length n which must contain the
!! diagonal elements of the matrix d.
real(wp), intent(in) :: Qtb(n) !! an input array of length n which must contain the first
!! n elements of the vector (q transpose)*b.
real(wp), intent(out) :: x(n) !! an output array of length n which contains the least
!! squares solution of the system a*x = b, sqrt(par)*d*x = 0,
!! for the output par.
real(wp), intent(out) :: Sdiag(n) !! an output array of length n which contains the
!! diagonal elements of the upper triangular matrix s.
real(wp), intent(inout) :: Wa1(n) !! work array of length n.
real(wp), intent(inout) :: Wa2(n) !! work array of length n.
integer :: i, iter, j, jm1, jp1, k, l, nsing
real(wp) :: dxnorm, fp, gnorm, parc, parl, paru, sum, temp
real(wp), parameter :: p1 = 1.0e-1_wp
real(wp), parameter :: p001 = 1.0e-3_wp
real(wp), parameter :: dwarf = dpmpar(2) !! the smallest positive magnitude
! 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) == zero .and. nsing == n) nsing = j - 1
if (nsing < n) Wa1(j) = zero
end do
if (nsing >= 1) then
do k = 1, nsing
j = nsing - k + 1
Wa1(j) = Wa1(j)/r(j, j)
temp = Wa1(j)
jm1 = j - 1
if (jm1 >= 1) then
do i = 1, jm1
Wa1(i) = Wa1(i) - r(i, j)*temp
end do
end if
end do
end if
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
do j = 1, n
Wa2(j) = Diag(j)*x(j)
end do
dxnorm = enorm(n, Wa2)
fp = dxnorm - Delta
if (fp <= p1*Delta) then
! termination.
if (iter == 0) Par = zero
else
! 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 = zero
if (nsing >= n) then
do j = 1, n
l = Ipvt(j)
Wa1(j) = Diag(l)*(Wa2(l)/dxnorm)
end do
do j = 1, n
sum = zero
jm1 = j - 1
if (jm1 >= 1) then
do i = 1, jm1
sum = sum + r(i, j)*Wa1(i)
end do
end if
Wa1(j) = (Wa1(j) - sum)/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
sum = zero
do i = 1, j
sum = sum + r(i, j)*Qtb(i)
end do
l = Ipvt(j)
Wa1(j) = sum/Diag(l)
end do
gnorm = enorm(n, Wa1)
paru = gnorm/Delta
if (paru == zero) paru = dwarf/min(Delta, p1)
! 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 == zero) Par = gnorm/dxnorm
! beginning of an iteration.
do
iter = iter + 1
! evaluate the function at the current value of par.
if (Par == zero) Par = max(dwarf, p001*paru)
temp = sqrt(Par)
do j = 1, n
Wa1(j) = temp*Diag(j)
end do
call qrsolv(n, r, Ldr, Ipvt, Wa1, Qtb, x, Sdiag, Wa2)
do j = 1, n
Wa2(j) = Diag(j)*x(j)
end do
dxnorm = enorm(n, Wa2)
temp = fp
fp = dxnorm - Delta
! if the function is small enough, accept the current value
! of par. also test for the exceptional cases where parl
! is zero or the number of iterations has reached 10.
if (abs(fp) <= p1*Delta .or. parl == zero .and. fp <= temp .and. &
temp < zero .or. iter == 10) then
if (iter == 0) Par = zero
exit
else
! 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)
jp1 = j + 1
if (n >= jp1) then
do i = jp1, n
Wa1(i) = Wa1(i) - r(i, j)*temp
end do
end if
end do
temp = enorm(n, Wa1)
parc = ((fp/Delta)/temp)/temp
! depending on the sign of the function, update parl or paru.
if (fp > zero) parl = max(parl, Par)
if (fp < zero) paru = min(paru, Par)
! compute an improved estimate for par.
Par = max(parl, Par + parc)
end if
end do ! end of an iteration.
end if
end subroutine lmpar