relibc_openlibm/slatec/qrsolv.f

*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