1*DECK QRSOLV
2      SUBROUTINE QRSOLV (N, R, LDR, IPVT, DIAG, QTB, X, SIGMA, WA)
3C***BEGIN PROLOGUE  QRSOLV
4C***SUBSIDIARY
5C***PURPOSE  Subsidiary to SNLS1 and SNLS1E
6C***LIBRARY   SLATEC
7C***TYPE      SINGLE PRECISION (QRSOLV-S, DQRSLV-D)
8C***AUTHOR  (UNKNOWN)
9C***DESCRIPTION
10C
11C     Given an M by N matrix A, an N by N diagonal matrix D,
12C     and an M-vector B, the problem is to determine an X which
13C     solves the system
14C
15C           A*X = B ,     D*X = 0 ,
16C
17C     in the least squares sense.
18C
19C     This subroutine completes the solution of the problem
20C     if it is provided with the necessary information from the
21C     QR factorization, with column pivoting, of A. That is, if
22C     A*P = Q*R, where P is a permutation matrix, Q has orthogonal
23C     columns, and R is an upper triangular matrix with diagonal
24C     elements of nonincreasing magnitude, then QRSOLV expects
25C     the full upper triangle of R, the permutation matrix P,
26C     and the first N components of (Q TRANSPOSE)*B. The system
27C     A*X = B, D*X = 0, is then equivalent to
28C
29C                  T       T
30C           R*Z = Q *B ,  P *D*P*Z = 0 ,
31C
32C     where X = P*Z. If this system does not have full rank,
33C     then a least squares solution is obtained. On output QRSOLV
34C     also provides an upper triangular matrix S such that
35C
36C            T   T               T
37C           P *(A *A + D*D)*P = S *S .
38C
39C     S is computed within QRSOLV and may be of separate interest.
40C
41C     The subroutine statement is
42C
43C       SUBROUTINE QRSOLV(N,R,LDR,IPVT,DIAG,QTB,X,SIGMA,WA)
44C
45C     where
46C
47C       N is a positive integer input variable set to the order of R.
48C
49C       R is an N by N array. On input the full upper triangle
50C         must contain the full upper triangle of the matrix R.
51C         On output the full upper triangle is unaltered, and the
52C         strict lower triangle contains the strict upper triangle
53C         (transposed) of the upper triangular matrix S.
54C
55C       LDR is a positive integer input variable not less than N
56C         which specifies the leading dimension of the array R.
57C
58C       IPVT is an integer input array of length N which defines the
59C         permutation matrix P such that A*P = Q*R. Column J of P
60C         is column IPVT(J) of the identity matrix.
61C
62C       DIAG is an input array of length N which must contain the
63C         diagonal elements of the matrix D.
64C
65C       QTB is an input array of length N which must contain the first
66C         N elements of the vector (Q TRANSPOSE)*B.
67C
68C       X is an output array of length N which contains the least
69C         squares solution of the system A*X = B, D*X = 0.
70C
71C       SIGMA is an output array of length N which contains the
72C         diagonal elements of the upper triangular matrix S.
73C
74C       WA is a work array of length N.
75C
76C***SEE ALSO  SNLS1, SNLS1E
77C***ROUTINES CALLED  (NONE)
78C***REVISION HISTORY  (YYMMDD)
79C   800301  DATE WRITTEN
80C   890831  Modified array declarations.  (WRB)
81C   891214  Prologue converted to Version 4.0 format.  (BAB)
82C   900326  Removed duplicate information from DESCRIPTION section.
83C           (WRB)
84C   900328  Added TYPE section.  (WRB)
85C***END PROLOGUE  QRSOLV
86      INTEGER N,LDR
87      INTEGER IPVT(*)
88      REAL R(LDR,*),DIAG(*),QTB(*),X(*),SIGMA(*),WA(*)
89      INTEGER I,J,JP1,K,KP1,L,NSING
90      REAL COS,COTAN,P5,P25,QTBPJ,SIN,SUM,TAN,TEMP,ZERO
91      SAVE P5, P25, ZERO
92      DATA P5,P25,ZERO /5.0E-1,2.5E-1,0.0E0/
93C***FIRST EXECUTABLE STATEMENT  QRSOLV
94      DO 20 J = 1, N
95         DO 10 I = J, N
96            R(I,J) = R(J,I)
97   10       CONTINUE
98         X(J) = R(J,J)
99         WA(J) = QTB(J)
100   20    CONTINUE
101C
102C     ELIMINATE THE DIAGONAL MATRIX D USING A GIVENS ROTATION.
103C
104      DO 100 J = 1, N
105C
106C        PREPARE THE ROW OF D TO BE ELIMINATED, LOCATING THE
107C        DIAGONAL ELEMENT USING P FROM THE QR FACTORIZATION.
108C
109         L = IPVT(J)
110         IF (DIAG(L) .EQ. ZERO) GO TO 90
111         DO 30 K = J, N
112            SIGMA(K) = ZERO
113   30       CONTINUE
114         SIGMA(J) = DIAG(L)
115C
116C        THE TRANSFORMATIONS TO ELIMINATE THE ROW OF D
117C        MODIFY ONLY A SINGLE ELEMENT OF (Q TRANSPOSE)*B
118C        BEYOND THE FIRST N, WHICH IS INITIALLY ZERO.
119C
120         QTBPJ = ZERO
121         DO 80 K = J, N
122C
123C           DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
124C           APPROPRIATE ELEMENT IN THE CURRENT ROW OF D.
125C
126            IF (SIGMA(K) .EQ. ZERO) GO TO 70
127            IF (ABS(R(K,K)) .GE. ABS(SIGMA(K))) GO TO 40
128               COTAN = R(K,K)/SIGMA(K)
129               SIN = P5/SQRT(P25+P25*COTAN**2)
130               COS = SIN*COTAN
131               GO TO 50
132   40       CONTINUE
133               TAN = SIGMA(K)/R(K,K)
134               COS = P5/SQRT(P25+P25*TAN**2)
135               SIN = COS*TAN
136   50       CONTINUE
137C
138C           COMPUTE THE MODIFIED DIAGONAL ELEMENT OF R AND
139C           THE MODIFIED ELEMENT OF ((Q TRANSPOSE)*B,0).
140C
141            R(K,K) = COS*R(K,K) + SIN*SIGMA(K)
142            TEMP = COS*WA(K) + SIN*QTBPJ
143            QTBPJ = -SIN*WA(K) + COS*QTBPJ
144            WA(K) = TEMP
145C
146C           ACCUMULATE THE TRANSFORMATION IN THE ROW OF S.
147C
148            KP1 = K + 1
149            IF (N .LT. KP1) GO TO 70
150            DO 60 I = KP1, N
151               TEMP = COS*R(I,K) + SIN*SIGMA(I)
152               SIGMA(I) = -SIN*R(I,K) + COS*SIGMA(I)
153               R(I,K) = TEMP
154   60          CONTINUE
155   70       CONTINUE
156   80       CONTINUE
157   90    CONTINUE
158C
159C        STORE THE DIAGONAL ELEMENT OF S AND RESTORE
160C        THE CORRESPONDING DIAGONAL ELEMENT OF R.
161C
162         SIGMA(J) = R(J,J)
163         R(J,J) = X(J)
164  100    CONTINUE
165C
166C     SOLVE THE TRIANGULAR SYSTEM FOR Z. IF THE SYSTEM IS
167C     SINGULAR, THEN OBTAIN A LEAST SQUARES SOLUTION.
168C
169      NSING = N
170      DO 110 J = 1, N
171         IF (SIGMA(J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1
172         IF (NSING .LT. N) WA(J) = ZERO
173  110    CONTINUE
174      IF (NSING .LT. 1) GO TO 150
175      DO 140 K = 1, NSING
176         J = NSING - K + 1
177         SUM = ZERO
178         JP1 = J + 1
179         IF (NSING .LT. JP1) GO TO 130
180         DO 120 I = JP1, NSING
181            SUM = SUM + R(I,J)*WA(I)
182  120       CONTINUE
183  130    CONTINUE
184         WA(J) = (WA(J) - SUM)/SIGMA(J)
185  140    CONTINUE
186  150 CONTINUE
187C
188C     PERMUTE THE COMPONENTS OF Z BACK TO COMPONENTS OF X.
189C
190      DO 160 J = 1, N
191         L = IPVT(J)
192         X(L) = WA(J)
193  160    CONTINUE
194      RETURN
195C
196C     LAST CARD OF SUBROUTINE QRSOLV.
197C
198      END
199