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