SUBROUTINE R1UPDT(m,n,s,ls,u,v,w,sing)
! 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.
!
! The subroutine statement is
!
! subroutine r1updt(m,n,s,ls,u,v,w,sing)
!
! where
!
! M is a positive integer input variable set to the number
! of rows of S.
!
! N is a positive integer input variable set to the number
! of columns of S. N must not exceed M.
!
! S is 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.
!
! LS is a positive integer input variable not less than
! (N*(2*M-N+1))/2.
!
! U is an input array of length M which must contain the
! vector U.
!
! V is 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.
!
! W is an output array of length M. W(I) contains information
! necessary to recover the Givens rotation GW(I) described
! above.
!
! SING is 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.
!
! Argonne National Laboratory. Minpack project. March 1980.
! Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More,
! John L. Nazareth
IMPLICIT NONE
INTEGER m,n,ls,i,j,jj,l,nmj,nm1
LOGICAL sing
real(dp) s(ls),u(m),v(n),w(m),cos1,cotan,giant,one
real(dp) p5,p25,sin1,tan1,tau,temp,zero
one = 1.0D0
p5 = 0.5D0
p25 = 0.25D0
zero = 0.0D0
! giant is the largest magnitude in the computer's arithmetic range
giant = spmpar(3)
! 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) goto 70
do nmj = 1, nm1
j = n - nmj
jj = jj - (m - j + 1)
w(j) = zero
if (v(j) == zero) goto 50
! Determine a givens rotation which eliminates the
! j-th element of v.
if (abs(v(n)) >= abs(v(j))) goto 20
cotan = v(n)/v(j)
sin1 = p5/sqrt(p25+p25*cotan**2)
cos1 = sin1*cotan
tau = one
if (abs(cos1)*giant > one) tau = one/cos1
goto 30
20 continue
tan1 = v(j)/v(n)
cos1 = p5/sqrt(p25+p25*tan1**2)
sin1 = cos1*tan1
tau = sin1
30 continue
! Apply the transformation to v and store the information
! necessary to recover the givens rotation.
v(n) = sin1*v(j) + cos1*v(n)
v(j) = tau
! Apply the transformation to s and extend the spike in w.
l = jj
do i = j, m
temp = cos1*s(l) - sin1*w(i)
w(i) = sin1*s(l) + cos1*w(i)
s(l) = temp
l = l + 1
end do
50 continue
end do
70 continue
! 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) goto 140
do j = 1, nm1
if (w(j) == zero) goto 120
! Determine a givens rotation which eliminates the
! j-th element of the spike.
if (abs(s(jj)) >= abs(w(j))) goto 90
cotan = s(jj)/w(j)
sin1 = p5/sqrt(p25+p25*cotan**2)
cos1 = sin1*cotan
tau = one
if ((abs(cos1)*giant) > one) tau = one/cos1
goto 100
90 continue
tan1 = w(j)/s(jj)
cos1 = p5/sqrt(p25+p25*tan1**2)
sin1 = cos1*tan1
tau = sin1
100 continue
! Apply the transformation to s and reduce the spike in w.
l = jj
do i = j, m
temp = cos1*s(l) + sin1*w(i)
w(i) = -sin1*s(l) + cos1*w(i)
s(l) = temp
l = l + 1
end do
! Store the information necessary to recover the
! givens rotation.
w(j) = tau
120 continue
! Test for zero diagonal elements in the output s.
if (s(jj) == zero) sing = .true.
jj = jj + (m - j + 1)
end do
140 continue
! 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.
return
end SUBROUTINE R1UPDT