+++ /dev/null
-*DECK QRSOLV
- SUBROUTINE QRSOLV (N, R, LDR, IPVT, DIAG, QTB, X, SIGMA, WA)
-C***BEGIN PROLOGUE QRSOLV
-C***SUBSIDIARY
-C***PURPOSE Subsidiary to SNLS1 and SNLS1E
-C***LIBRARY SLATEC
-C***TYPE SINGLE PRECISION (QRSOLV-S, DQRSLV-D)
-C***AUTHOR (UNKNOWN)
-C***DESCRIPTION
-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 QRSOLV 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 QRSOLV
-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 QRSOLV and may be of separate interest.
-C
-C The subroutine statement is
-C
-C SUBROUTINE QRSOLV(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 SNLS1, SNLS1E
-C***ROUTINES CALLED (NONE)
-C***REVISION HISTORY (YYMMDD)
-C 800301 DATE WRITTEN
-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 QRSOLV
- INTEGER N,LDR
- INTEGER IPVT(*)
- REAL R(LDR,*),DIAG(*),QTB(*),X(*),SIGMA(*),WA(*)
- INTEGER I,J,JP1,K,KP1,L,NSING
- REAL COS,COTAN,P5,P25,QTBPJ,SIN,SUM,TAN,TEMP,ZERO
- SAVE P5, P25, ZERO
- DATA P5,P25,ZERO /5.0E-1,2.5E-1,0.0E0/
-C***FIRST EXECUTABLE STATEMENT QRSOLV
- 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 QRSOLV.
-C
- END
git://source.jalview.org / jabaws.git / blobdiff