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(2m-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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | m | ||||
| integer | :: | n | ||||
| real(kind=wp) | :: | s(ls) | ||||
| integer | :: | ls | ||||
| real(kind=wp) | :: | u(m) | ||||
| real(kind=wp) | :: | v(n) | ||||
| real(kind=wp) | :: | w(m) | ||||
| logical | :: | sing |
subroutine r1updt(m,n,s,ls,u,v,w,sing) implicit none integer m , n , ls logical sing real(wp) s(ls) , u(m) , v(n) , w(m) integer i , j , jj , l , nmj , nm1 real(wp) cos , cotan , giant , sin , & & tan , tau , temp real(wp),parameter :: p5 = 5.0e-1_wp real(wp),parameter :: p25 = 2.5e-1_wp ! ! giant is the largest magnitude. ! giant = dpmpar(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 enddo ! ! 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 endif ! ! 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 enddo endif enddo endif ! ! add the spike from the rank 1 update to w. ! do i = 1 , m w(i) = w(i) + v(n)*u(i) enddo ! ! 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 endif ! ! 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 enddo ! ! store the information necessary to recover the ! givens rotation. ! w(j) = tau endif ! ! test for zero diagonal elements in the output s. ! if ( s(jj)==zero ) sing = .true. jj = jj + (m-j+1) enddo endif ! ! move w back into the last column of the output s. ! l = jj do i = n , m s(l) = w(i) l = l + 1 enddo if ( s(jj)==zero ) sing = .true. end subroutine r1updt