2 SUBROUTINE RWUPDT (N, R, LDR, W, B, ALPHA, COS, SIN)
3 C***BEGIN PROLOGUE RWUPDT
5 C***PURPOSE Subsidiary to SNLS1 and SNLS1E
7 C***TYPE SINGLE PRECISION (RWUPDT-S, DWUPDT-D)
11 C Given an N by N upper triangular matrix R, this subroutine
12 C computes the QR decomposition of the matrix formed when a row
13 C is added to R. If the row is specified by the vector W, then
14 C RWUPDT determines an orthogonal matrix Q such that when the
15 C N+1 by N matrix composed of R augmented by W is premultiplied
16 C by (Q TRANSPOSE), the resulting matrix is upper trapezoidal.
17 C The orthogonal matrix Q is the product of N transformations
19 C G(1)*G(2)* ... *G(N)
21 C where G(I) is a Givens rotation in the (I,N+1) plane which
22 C eliminates elements in the I-th plane. RWUPDT also
23 C computes the product (Q TRANSPOSE)*C where C is the
24 C (N+1)-vector (b,alpha). Q itself is not accumulated, rather
25 C the information to recover the G rotations is supplied.
27 C The subroutine statement is
29 C SUBROUTINE RWUPDT(N,R,LDR,W,B,ALPHA,COS,SIN)
33 C N is a positive integer input variable set to the order of R.
35 C R is an N by N array. On input the upper triangular part of
36 C R must contain the matrix to be updated. On output R
37 C contains the updated triangular matrix.
39 C LDR is a positive integer input variable not less than N
40 C which specifies the leading dimension of the array R.
42 C W is an input array of length N which must contain the row
43 C vector to be added to R.
45 C B is an array of length N. On input B must contain the
46 C first N elements of the vector C. On output B contains
47 C the first N elements of the vector (Q TRANSPOSE)*C.
49 C ALPHA is a variable. On input ALPHA must contain the
50 C (N+1)-st element of the vector C. On output ALPHA contains
51 C the (N+1)-st element of the vector (Q TRANSPOSE)*C.
53 C COS is an output array of length N which contains the
54 C cosines of the transforming Givens rotations.
56 C SIN is an output array of length N which contains the
57 C sines of the transforming Givens rotations.
59 C***SEE ALSO SNLS1, SNLS1E
60 C***ROUTINES CALLED (NONE)
61 C***REVISION HISTORY (YYMMDD)
63 C 890831 Modified array declarations. (WRB)
64 C 891214 Prologue converted to Version 4.0 format. (BAB)
65 C 900326 Removed duplicate information from DESCRIPTION section.
67 C 900328 Added TYPE section. (WRB)
68 C***END PROLOGUE RWUPDT
71 REAL R(LDR,*),W(*),B(*),COS(*),SIN(*)
73 REAL COTAN,ONE,P5,P25,ROWJ,TAN,TEMP,ZERO
74 SAVE ONE, P5, P25, ZERO
75 DATA ONE,P5,P25,ZERO /1.0E0,5.0E-1,2.5E-1,0.0E0/
76 C***FIRST EXECUTABLE STATEMENT RWUPDT
81 C APPLY THE PREVIOUS TRANSFORMATIONS TO
82 C R(I,J), I=1,2,...,J-1, AND TO W(J).
84 IF (JM1 .LT. 1) GO TO 20
86 TEMP = COS(I)*R(I,J) + SIN(I)*ROWJ
87 ROWJ = -SIN(I)*R(I,J) + COS(I)*ROWJ
92 C DETERMINE A GIVENS ROTATION WHICH ELIMINATES W(J).
96 IF (ROWJ .EQ. ZERO) GO TO 50
97 IF (ABS(R(J,J)) .GE. ABS(ROWJ)) GO TO 30
99 SIN(J) = P5/SQRT(P25+P25*COTAN**2)
100 COS(J) = SIN(J)*COTAN
104 COS(J) = P5/SQRT(P25+P25*TAN**2)
108 C APPLY THE CURRENT TRANSFORMATION TO R(J,J), B(J), AND ALPHA.
110 R(J,J) = COS(J)*R(J,J) + SIN(J)*ROWJ
111 TEMP = COS(J)*B(J) + SIN(J)*ALPHA
112 ALPHA = -SIN(J)*B(J) + COS(J)*ALPHA
118 C LAST CARD OF SUBROUTINE RWUPDT.