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
t
(s + u*v )*q
is again lower trapezoidal.
this subroutine 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, gw rotations is returned.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | m | a positive integer input variable set to the number of rows of s. |
||
| integer, | intent(in) | :: | n | a positive integer input variable set to the number of columns of s. n must not exceed m. |
||
| real(kind=wp), | intent(inout) | :: | s(Ls) | an array of length ls. on input s must contain the lower trapezoidal matrix s stored by columns. on output s contains the lower trapezoidal matrix produced as described above. |
||
| integer, | intent(in) | :: | Ls | a positive integer input variable not less than (n(2m-n+1))/2. |
||
| real(kind=wp), | intent(in) | :: | u(m) | an input array of length m which must contain the vector u. |
||
| real(kind=wp), | intent(inout) | :: | v(n) | an array of length n. on input v must contain the vector v. on output v(i) contains the information necessary to recover the givens rotation gv(i) described above. |
||
| real(kind=wp), | intent(out) | :: | w(m) | an output array of length m. w(i) contains information necessary to recover the givens rotation gw(i) described above. |
||
| logical, | intent(out) | :: | Sing | a logical output variable. sing is set true if any of the diagonal elements of the output s are zero. otherwise sing is set false. |
subroutine r1updt(m, n, s, Ls, u, v, w, Sing)
implicit none
integer, intent(in) :: m !! a positive integer input variable set to the number
!! of rows of s.
integer, intent(in) :: n !! a positive integer input variable set to the number
!! of columns of s. n must not exceed m.
integer, intent(in) :: Ls !! a positive integer input variable not less than
!! (n*(2*m-n+1))/2.
logical, intent(out) :: Sing !! a logical output variable. sing is set true if any
!! of the diagonal elements of the output s are zero. otherwise
!! sing is set false.
real(wp), intent(inout) :: s(Ls) !! an array of length ls. on input s must contain the lower
!! trapezoidal matrix s stored by columns. on output s contains
!! the lower trapezoidal matrix produced as described above.
real(wp), intent(in) :: u(m) !! an input array of length m which must contain the
!! vector u.
real(wp), intent(inout) :: v(n) !! an array of length n. on input v must contain the vector
!! v. on output v(i) contains the information necessary to
!! recover the givens rotation gv(i) described above.
real(wp), intent(out) :: w(m) !! an output array of length m. w(i) contains information
!! necessary to recover the givens rotation gw(i) described
!! above.
integer :: i, j, jj, l, nmj, nm1
real(wp) :: cos, cotan, sin, tan, tau, temp
real(wp), parameter :: p5 = 5.0e-1_wp
real(wp), parameter :: p25 = 2.5e-1_wp
real(wp), parameter :: giant = dpmpar(3) !! the largest magnitude.
! 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.
nm1 = n - 1
if (nm1 >= 1) then
do nmj = 1, nm1
j = n - nmj
jj = jj - (m - j + 1)
w(j) = zero
if (v(j) /= zero) then
! determine a givens rotation which eliminates the
! j-th element of v.
if (abs(v(n)) >= abs(v(j))) then
tan = v(j)/v(n)
cos = p5/sqrt(p25 + p25*tan**2)
sin = cos*tan
tau = sin
else
cotan = v(n)/v(j)
sin = p5/sqrt(p25 + p25*cotan**2)
cos = sin*cotan
tau = one
if (abs(cos)*giant > one) tau = one/cos
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
end if
! add the spike from the rank 1 update to w.
do i = 1, m
w(i) = w(i) + v(n)*u(i)
end do
! eliminate the spike.
Sing = .false.
if (nm1 >= 1) then
do j = 1, nm1
if (w(j) /= zero) then
! determine a givens rotation which eliminates the
! j-th element of the spike.
if (abs(s(jj)) >= abs(w(j))) then
tan = w(j)/s(jj)
cos = p5/sqrt(p25 + p25*tan**2)
sin = cos*tan
tau = sin
else
cotan = s(jj)/w(j)
sin = p5/sqrt(p25 + p25*cotan**2)
cos = sin*cotan
tau = one
if (abs(cos)*giant > one) tau = one/cos
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) == zero) Sing = .true.
jj = jj + (m - j + 1)
end do
end if
! 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) == zero) Sing = .true.
end subroutine r1updt