*DECK D1UPDT
SUBROUTINE D1UPDT (M, N, S, LS, U, V, W, SING)
C***BEGIN PROLOGUE D1UPDT
C***SUBSIDIARY
C***PURPOSE Subsidiary to DNSQ and DNSQE
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (R1UPDT-S, D1UPDT-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
C
C Given an M by N lower trapezoidal matrix S, an M-vector U,
C and an N-vector V, the problem is to determine an
C orthogonal matrix Q such that
C
C t
C (S + U*V )*Q
C
C is again lower trapezoidal.
C
C This subroutine determines Q as the product of 2*(N - 1)
C transformations
C
C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
C
C where GV(I), GW(I) are Givens rotations in the (I,N) plane
C which eliminate elements in the I-th and N-th planes,
C respectively. Q itself is not accumulated, rather the
C information to recover the GV, GW rotations is returned.
C
C The SUBROUTINE statement is
C
C SUBROUTINE D1UPDT(M,N,S,LS,U,V,W,SING)
C
C where
C
C M is a positive integer input variable set to the number
C of rows of S.
C
C N is a positive integer input variable set to the number
C of columns of S. N must not exceed M.
C
C S is an array of length LS. On input S must contain the lower
C trapezoidal matrix S stored by columns. On output S contains
C the lower trapezoidal matrix produced as described above.
C
C LS is a positive integer input variable not less than
C (N*(2*M-N+1))/2.
C
C U is an input array of length M which must contain the
C vector U.
C
C V is an array of length N. On input V must contain the vector
C V. On output V(I) contains the information necessary to
C recover the Givens rotation GV(I) described above.
C
C W is an output array of length M. W(I) contains information
C necessary to recover the Givens rotation GW(I) described
C above.
C
C SING is a LOGICAL output variable. SING is set TRUE if any
C of the diagonal elements of the output S are zero. Otherwise
C SING is set FALSE.
C
C***SEE ALSO DNSQ, DNSQE
C***ROUTINES CALLED D1MACH
C***REVISION HISTORY (YYMMDD)
C 800301 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 900328 Added TYPE section. (WRB)
C***END PROLOGUE D1UPDT
DOUBLE PRECISION D1MACH
INTEGER I, J, JJ, L, LS, M, N, NM1, NMJ
DOUBLE PRECISION COS, COTAN, GIANT, ONE, P25, P5, S(*),
1 SIN, TAN, TAU, TEMP, U(*), V(*), W(*), ZERO
LOGICAL SING
SAVE ONE, P5, P25, ZERO
DATA ONE,P5,P25,ZERO /1.0D0,5.0D-1,2.5D-1,0.0D0/
C
C GIANT IS THE LARGEST MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT D1UPDT
GIANT = D1MACH(2)
C
C INITIALIZE THE DIAGONAL ELEMENT POINTER.
C
JJ = (N*(2*M - N + 1))/2 - (M - N)
C
C MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W.
C
L = JJ
DO 10 I = N, M
W(I) = S(L)
L = L + 1
10 CONTINUE
C
C ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR
C IN SUCH A WAY THAT A SPIKE IS INTRODUCED INTO W.
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 70
DO 60 NMJ = 1, NM1
J = N - NMJ
JJ = JJ - (M - J + 1)
W(J) = ZERO
IF (V(J) .EQ. ZERO) GO TO 50
C
C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
C J-TH ELEMENT OF V.
C
IF (ABS(V(N)) .GE. ABS(V(J))) GO TO 20
COTAN = V(N)/V(J)
SIN = P5/SQRT(P25+P25*COTAN**2)
COS = SIN*COTAN
TAU = ONE
IF (ABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
GO TO 30
20 CONTINUE
TAN = V(J)/V(N)
COS = P5/SQRT(P25+P25*TAN**2)
SIN = COS*TAN
TAU = SIN
30 CONTINUE
C
C APPLY THE TRANSFORMATION TO V AND STORE THE INFORMATION
C NECESSARY TO RECOVER THE GIVENS ROTATION.
C
V(N) = SIN*V(J) + COS*V(N)
V(J) = TAU
C
C APPLY THE TRANSFORMATION TO S AND EXTEND THE SPIKE IN W.
C
L = JJ
DO 40 I = J, M
TEMP = COS*S(L) - SIN*W(I)
W(I) = SIN*S(L) + COS*W(I)
S(L) = TEMP
L = L + 1
40 CONTINUE
50 CONTINUE
60 CONTINUE
70 CONTINUE
C
C ADD THE SPIKE FROM THE RANK 1 UPDATE TO W.
C
DO 80 I = 1, M
W(I) = W(I) + V(N)*U(I)
80 CONTINUE
C
C ELIMINATE THE SPIKE.
C
SING = .FALSE.
IF (NM1 .LT. 1) GO TO 140
DO 130 J = 1, NM1
IF (W(J) .EQ. ZERO) GO TO 120
C
C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
C J-TH ELEMENT OF THE SPIKE.
C
IF (ABS(S(JJ)) .GE. ABS(W(J))) GO TO 90
COTAN = S(JJ)/W(J)
SIN = P5/SQRT(P25+P25*COTAN**2)
COS = SIN*COTAN
TAU = ONE
IF (ABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
GO TO 100
90 CONTINUE
TAN = W(J)/S(JJ)
COS = P5/SQRT(P25+P25*TAN**2)
SIN = COS*TAN
TAU = SIN
100 CONTINUE
C
C APPLY THE TRANSFORMATION TO S AND REDUCE THE SPIKE IN W.
C
L = JJ
DO 110 I = J, M
TEMP = COS*S(L) + SIN*W(I)
W(I) = -SIN*S(L) + COS*W(I)
S(L) = TEMP
L = L + 1
110 CONTINUE
C
C STORE THE INFORMATION NECESSARY TO RECOVER THE
C GIVENS ROTATION.
C
W(J) = TAU
120 CONTINUE
C
C TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S.
C
IF (S(JJ) .EQ. ZERO) SING = .TRUE.
JJ = JJ + (M - J + 1)
130 CONTINUE
140 CONTINUE
C
C MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S.
C
L = JJ
DO 150 I = N, M
S(L) = W(I)
L = L + 1
150 CONTINUE
IF (S(JJ) .EQ. ZERO) SING = .TRUE.
RETURN
C
C LAST CARD OF SUBROUTINE D1UPDT.
C
END