*DECK DQRSLV
SUBROUTINE DQRSLV (N, R, LDR, IPVT, DIAG, QTB, X, SIGMA, WA)
C***BEGIN PROLOGUE DQRSLV
C***SUBSIDIARY
C***PURPOSE Subsidiary to DNLS1 and DNLS1E
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (QRSOLV-S, DQRSLV-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
C
C **** Double Precision version of QRSOLV ****
C
C Given an M by N matrix A, an N by N diagonal matrix D,
C and an M-vector B, the problem is to determine an X which
C solves the system
C
C A*X = B , D*X = 0 ,
C
C in the least squares sense.
C
C This subroutine completes the solution of the problem
C if it is provided with the necessary information from the
C QR factorization, with column pivoting, of A. That is, if
C A*P = Q*R, where P is a permutation matrix, Q has orthogonal
C columns, and R is an upper triangular matrix with diagonal
C elements of nonincreasing magnitude, then DQRSLV expects
C the full upper triangle of R, the permutation matrix P,
C and the first N components of (Q TRANSPOSE)*B. The system
C A*X = B, D*X = 0, is then equivalent to
C
C T T
C R*Z = Q *B , P *D*P*Z = 0 ,
C
C where X = P*Z. If this system does not have full rank,
C then a least squares solution is obtained. On output DQRSLV
C also provides an upper triangular matrix S such that
C
C T T T
C P *(A *A + D*D)*P = S *S .
C
C S is computed within DQRSLV and may be of separate interest.
C
C The subroutine statement is
C
C SUBROUTINE DQRSLV(N,R,LDR,IPVT,DIAG,QTB,X,SIGMA,WA)
C
C where
C
C N is a positive integer input variable set to the order of R.
C
C R is an N by N array. On input the full upper triangle
C must contain the full upper triangle of the matrix R.
C On output the full upper triangle is unaltered, and the
C strict lower triangle contains the strict upper triangle
C (transposed) of the upper triangular matrix S.
C
C LDR is a positive integer input variable not less than N
C which specifies the leading dimension of the array R.
C
C IPVT is an integer input array of length N which defines the
C permutation matrix P such that A*P = Q*R. Column J of P
C is column IPVT(J) of the identity matrix.
C
C DIAG is an input array of length N which must contain the
C diagonal elements of the matrix D.
C
C QTB is an input array of length N which must contain the first
C N elements of the vector (Q TRANSPOSE)*B.
C
C X is an output array of length N which contains the least
C squares solution of the system A*X = B, D*X = 0.
C
C SIGMA is an output array of length N which contains the
C diagonal elements of the upper triangular matrix S.
C
C WA is a work array of length N.
C
C***SEE ALSO DNLS1, DNLS1E
C***ROUTINES CALLED (NONE)
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 DQRSLV
INTEGER N,LDR
INTEGER IPVT(*)
DOUBLE PRECISION R(LDR,*),DIAG(*),QTB(*),X(*),SIGMA(*),WA(*)
INTEGER I,J,JP1,K,KP1,L,NSING
DOUBLE PRECISION COS,COTAN,P5,P25,QTBPJ,SIN,SUM,TAN,TEMP,ZERO
SAVE P5, P25, ZERO
DATA P5,P25,ZERO /5.0D-1,2.5D-1,0.0D0/
C***FIRST EXECUTABLE STATEMENT DQRSLV
DO 20 J = 1, N
DO 10 I = J, N
R(I,J) = R(J,I)
10 CONTINUE
X(J) = R(J,J)
WA(J) = QTB(J)
20 CONTINUE
C
C ELIMINATE THE DIAGONAL MATRIX D USING A GIVENS ROTATION.
C
DO 100 J = 1, N
C
C PREPARE THE ROW OF D TO BE ELIMINATED, LOCATING THE
C DIAGONAL ELEMENT USING P FROM THE QR FACTORIZATION.
C
L = IPVT(J)
IF (DIAG(L) .EQ. ZERO) GO TO 90
DO 30 K = J, N
SIGMA(K) = ZERO
30 CONTINUE
SIGMA(J) = DIAG(L)
C
C THE TRANSFORMATIONS TO ELIMINATE THE ROW OF D
C MODIFY ONLY A SINGLE ELEMENT OF (Q TRANSPOSE)*B
C BEYOND THE FIRST N, WHICH IS INITIALLY ZERO.
C
QTBPJ = ZERO
DO 80 K = J, N
C
C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
C APPROPRIATE ELEMENT IN THE CURRENT ROW OF D.
C
IF (SIGMA(K) .EQ. ZERO) GO TO 70
IF (ABS(R(K,K)) .GE. ABS(SIGMA(K))) GO TO 40
COTAN = R(K,K)/SIGMA(K)
SIN = P5/SQRT(P25+P25*COTAN**2)
COS = SIN*COTAN
GO TO 50
40 CONTINUE
TAN = SIGMA(K)/R(K,K)
COS = P5/SQRT(P25+P25*TAN**2)
SIN = COS*TAN
50 CONTINUE
C
C COMPUTE THE MODIFIED DIAGONAL ELEMENT OF R AND
C THE MODIFIED ELEMENT OF ((Q TRANSPOSE)*B,0).
C
R(K,K) = COS*R(K,K) + SIN*SIGMA(K)
TEMP = COS*WA(K) + SIN*QTBPJ
QTBPJ = -SIN*WA(K) + COS*QTBPJ
WA(K) = TEMP
C
C ACCUMULATE THE TRANSFORMATION IN THE ROW OF S.
C
KP1 = K + 1
IF (N .LT. KP1) GO TO 70
DO 60 I = KP1, N
TEMP = COS*R(I,K) + SIN*SIGMA(I)
SIGMA(I) = -SIN*R(I,K) + COS*SIGMA(I)
R(I,K) = TEMP
60 CONTINUE
70 CONTINUE
80 CONTINUE
90 CONTINUE
C
C STORE THE DIAGONAL ELEMENT OF S AND RESTORE
C THE CORRESPONDING DIAGONAL ELEMENT OF R.
C
SIGMA(J) = R(J,J)
R(J,J) = X(J)
100 CONTINUE
C
C SOLVE THE TRIANGULAR SYSTEM FOR Z. IF THE SYSTEM IS
C SINGULAR, THEN OBTAIN A LEAST SQUARES SOLUTION.
C
NSING = N
DO 110 J = 1, N
IF (SIGMA(J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1
IF (NSING .LT. N) WA(J) = ZERO
110 CONTINUE
IF (NSING .LT. 1) GO TO 150
DO 140 K = 1, NSING
J = NSING - K + 1
SUM = ZERO
JP1 = J + 1
IF (NSING .LT. JP1) GO TO 130
DO 120 I = JP1, NSING
SUM = SUM + R(I,J)*WA(I)
120 CONTINUE
130 CONTINUE
WA(J) = (WA(J) - SUM)/SIGMA(J)
140 CONTINUE
150 CONTINUE
C
C PERMUTE THE COMPONENTS OF Z BACK TO COMPONENTS OF X.
C
DO 160 J = 1, N
L = IPVT(J)
X(L) = WA(J)
160 CONTINUE
RETURN
C
C LAST CARD OF SUBROUTINE DQRSLV.
C
END
