binaries/src/disembl/Tisean_3.0.1/source_f/slatec/qrsolv.f

   1 *DECK QRSOLV
   2       SUBROUTINE QRSOLV (N, R, LDR, IPVT, DIAG, QTB, X, SIGMA, WA)
   3 C***BEGIN PROLOGUE  QRSOLV
   4 C***SUBSIDIARY
   5 C***PURPOSE  Subsidiary to SNLS1 and SNLS1E
   6 C***LIBRARY   SLATEC
   7 C***TYPE      SINGLE PRECISION (QRSOLV-S, DQRSLV-D)
   8 C***AUTHOR  (UNKNOWN)
   9 C***DESCRIPTION
  10 C
  11 C     Given an M by N matrix A, an N by N diagonal matrix D,
  12 C     and an M-vector B, the problem is to determine an X which
  13 C     solves the system
  14 C
  15 C           A*X = B ,     D*X = 0 ,
  16 C
  17 C     in the least squares sense.
  18 C
  19 C     This subroutine completes the solution of the problem
  20 C     if it is provided with the necessary information from the
  21 C     QR factorization, with column pivoting, of A. That is, if
  22 C     A*P = Q*R, where P is a permutation matrix, Q has orthogonal
  23 C     columns, and R is an upper triangular matrix with diagonal
  24 C     elements of nonincreasing magnitude, then QRSOLV expects
  25 C     the full upper triangle of R, the permutation matrix P,
  26 C     and the first N components of (Q TRANSPOSE)*B. The system
  27 C     A*X = B, D*X = 0, is then equivalent to
  28 C
  29 C                  T       T
  30 C           R*Z = Q *B ,  P *D*P*Z = 0 ,
  31 C
  32 C     where X = P*Z. If this system does not have full rank,
  33 C     then a least squares solution is obtained. On output QRSOLV
  34 C     also provides an upper triangular matrix S such that
  35 C
  36 C            T   T               T
  37 C           P *(A *A + D*D)*P = S *S .
  38 C
  39 C     S is computed within QRSOLV and may be of separate interest.
  40 C
  41 C     The subroutine statement is
  42 C
  43 C       SUBROUTINE QRSOLV(N,R,LDR,IPVT,DIAG,QTB,X,SIGMA,WA)
  44 C
  45 C     where
  46 C
  47 C       N is a positive integer input variable set to the order of R.
  48 C
  49 C       R is an N by N array. On input the full upper triangle
  50 C         must contain the full upper triangle of the matrix R.
  51 C         On output the full upper triangle is unaltered, and the
  52 C         strict lower triangle contains the strict upper triangle
  53 C         (transposed) of the upper triangular matrix S.
  54 C
  55 C       LDR is a positive integer input variable not less than N
  56 C         which specifies the leading dimension of the array R.
  57 C
  58 C       IPVT is an integer input array of length N which defines the
  59 C         permutation matrix P such that A*P = Q*R. Column J of P
  60 C         is column IPVT(J) of the identity matrix.
  61 C
  62 C       DIAG is an input array of length N which must contain the
  63 C         diagonal elements of the matrix D.
  64 C
  65 C       QTB is an input array of length N which must contain the first
  66 C         N elements of the vector (Q TRANSPOSE)*B.
  67 C
  68 C       X is an output array of length N which contains the least
  69 C         squares solution of the system A*X = B, D*X = 0.
  70 C
  71 C       SIGMA is an output array of length N which contains the
  72 C         diagonal elements of the upper triangular matrix S.
  73 C
  74 C       WA is a work array of length N.
  75 C
  76 C***SEE ALSO  SNLS1, SNLS1E
  77 C***ROUTINES CALLED  (NONE)
  78 C***REVISION HISTORY  (YYMMDD)
  79 C   800301  DATE WRITTEN
  80 C   890831  Modified array declarations.  (WRB)
  81 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  82 C   900326  Removed duplicate information from DESCRIPTION section.
  83 C           (WRB)
  84 C   900328  Added TYPE section.  (WRB)
  85 C***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/
  93 C***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
 101 C
 102 C     ELIMINATE THE DIAGONAL MATRIX D USING A GIVENS ROTATION.
 103 C
 104       DO 100 J = 1, N
 105 C
 106 C        PREPARE THE ROW OF D TO BE ELIMINATED, LOCATING THE
 107 C        DIAGONAL ELEMENT USING P FROM THE QR FACTORIZATION.
 108 C
 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)
 115 C
 116 C        THE TRANSFORMATIONS TO ELIMINATE THE ROW OF D
 117 C        MODIFY ONLY A SINGLE ELEMENT OF (Q TRANSPOSE)*B
 118 C        BEYOND THE FIRST N, WHICH IS INITIALLY ZERO.
 119 C
 120          QTBPJ = ZERO
 121          DO 80 K = J, N
 122 C
 123 C           DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
 124 C           APPROPRIATE ELEMENT IN THE CURRENT ROW OF D.
 125 C
 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
 137 C
 138 C           COMPUTE THE MODIFIED DIAGONAL ELEMENT OF R AND
 139 C           THE MODIFIED ELEMENT OF ((Q TRANSPOSE)*B,0).
 140 C
 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
 145 C
 146 C           ACCUMULATE THE TRANSFORMATION IN THE ROW OF S.
 147 C
 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
 158 C
 159 C        STORE THE DIAGONAL ELEMENT OF S AND RESTORE
 160 C        THE CORRESPONDING DIAGONAL ELEMENT OF R.
 161 C
 162          SIGMA(J) = R(J,J)
 163          R(J,J) = X(J)
 164   100    CONTINUE
 165 C
 166 C     SOLVE THE TRIANGULAR SYSTEM FOR Z. IF THE SYSTEM IS
 167 C     SINGULAR, THEN OBTAIN A LEAST SQUARES SOLUTION.
 168 C
 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
 187 C
 188 C     PERMUTE THE COMPONENTS OF Z BACK TO COMPONENTS OF X.
 189 C
 190       DO 160 J = 1, N
 191          L = IPVT(J)
 192          X(L) = WA(J)
 193   160    CONTINUE
 194       RETURN
 195 C
 196 C     LAST CARD OF SUBROUTINE QRSOLV.
 197 C
 198       END