website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/qrfac.f

   1 *DECK QRFAC
   2       SUBROUTINE QRFAC (M, N, A, LDA, PIVOT, IPVT, LIPVT, SIGMA, ACNORM,
   3      +   WA)
   4 C***BEGIN PROLOGUE  QRFAC
   5 C***SUBSIDIARY
   6 C***PURPOSE  Subsidiary to SNLS1, SNLS1E, SNSQ and SNSQE
   7 C***LIBRARY   SLATEC
   8 C***TYPE      SINGLE PRECISION (QRFAC-S, DQRFAC-D)
   9 C***AUTHOR  (UNKNOWN)
  10 C***DESCRIPTION
  11 C
  12 C     This subroutine uses Householder transformations with column
  13 C     pivoting (optional) to compute a QR factorization of the
  14 C     M by N matrix A. That is, QRFAC determines an orthogonal
  15 C     matrix Q, a permutation matrix P, and an upper trapezoidal
  16 C     matrix R with diagonal elements of nonincreasing magnitude,
  17 C     such that A*P = Q*R. The Householder transformation for
  18 C     column K, K = 1,2,...,MIN(M,N), is of the form
  19 C
  20 C                           T
  21 C           I - (1/U(K))*U*U
  22 C
  23 C     where U has zeros in the first K-1 positions. The form of
  24 C     this transformation and the method of pivoting first
  25 C     appeared in the corresponding LINPACK subroutine.
  26 C
  27 C     The subroutine statement is
  28 C
  29 C       SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
  30 C
  31 C     where
  32 C
  33 C       M is a positive integer input variable set to the number
  34 C         of rows of A.
  35 C
  36 C       N is a positive integer input variable set to the number
  37 C         of columns of A.
  38 C
  39 C       A is an M by N array. On input A contains the matrix for
  40 C         which the QR factorization is to be computed. On output
  41 C         the strict upper trapezoidal part of A contains the strict
  42 C         upper trapezoidal part of R, and the lower trapezoidal
  43 C         part of A contains a factored form of Q (the non-trivial
  44 C         elements of the U vectors described above).
  45 C
  46 C       LDA is a positive integer input variable not less than M
  47 C         which specifies the leading dimension of the array A.
  48 C
  49 C       PIVOT is a logical input variable. If pivot is set .TRUE.,
  50 C         then column pivoting is enforced. If pivot is set .FALSE.,
  51 C         then no column pivoting is done.
  52 C
  53 C       IPVT is an integer output array of length LIPVT. IPVT
  54 C         defines the permutation matrix P such that A*P = Q*R.
  55 C         Column J of P is column IPVT(J) of the identity matrix.
  56 C         If pivot is .FALSE., IPVT is not referenced.
  57 C
  58 C       LIPVT is a positive integer input variable. If PIVOT is
  59 C             .FALSE., then LIPVT may be as small as 1. If PIVOT is
  60 C             .TRUE., then LIPVT must be at least N.
  61 C
  62 C       SIGMA is an output array of length N which contains the
  63 C         diagonal elements of R.
  64 C
  65 C       ACNORM is an output array of length N which contains the
  66 C         norms of the corresponding columns of the input matrix A.
  67 C         If this information is not needed, then ACNORM can coincide
  68 C         with SIGMA.
  69 C
  70 C       WA is a work array of length N. If pivot is .FALSE., then WA
  71 C         can coincide with SIGMA.
  72 C
  73 C***SEE ALSO  SNLS1, SNLS1E, SNSQ, SNSQE
  74 C***ROUTINES CALLED  ENORM, R1MACH
  75 C***REVISION HISTORY  (YYMMDD)
  76 C   800301  DATE WRITTEN
  77 C   890531  Changed all specific intrinsics to generic.  (WRB)
  78 C   890831  Modified array declarations.  (WRB)
  79 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  80 C   900326  Removed duplicate information from DESCRIPTION section.
  81 C           (WRB)
  82 C   900328  Added TYPE section.  (WRB)
  83 C***END PROLOGUE  QRFAC
  84       INTEGER M,N,LDA,LIPVT
  85       INTEGER IPVT(*)
  86       LOGICAL PIVOT
  87       REAL A(LDA,*),SIGMA(*),ACNORM(*),WA(*)
  88       INTEGER I,J,JP1,K,KMAX,MINMN
  89       REAL AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
  90       REAL R1MACH,ENORM
  91       SAVE ONE, P05, ZERO
  92       DATA ONE,P05,ZERO /1.0E0,5.0E-2,0.0E0/
  93 C***FIRST EXECUTABLE STATEMENT  QRFAC
  94       EPSMCH = R1MACH(4)
  95 C
  96 C     COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
  97 C
  98       DO 10 J = 1, N
  99          ACNORM(J) = ENORM(M,A(1,J))
 100          SIGMA(J) = ACNORM(J)
 101          WA(J) = SIGMA(J)
 102          IF (PIVOT) IPVT(J) = J
 103    10    CONTINUE
 104 C
 105 C     REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.
 106 C
 107       MINMN = MIN(M,N)
 108       DO 110 J = 1, MINMN
 109          IF (.NOT.PIVOT) GO TO 40
 110 C
 111 C        BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.
 112 C
 113          KMAX = J
 114          DO 20 K = J, N
 115             IF (SIGMA(K) .GT. SIGMA(KMAX)) KMAX = K
 116    20       CONTINUE
 117          IF (KMAX .EQ. J) GO TO 40
 118          DO 30 I = 1, M
 119             TEMP = A(I,J)
 120             A(I,J) = A(I,KMAX)
 121             A(I,KMAX) = TEMP
 122    30       CONTINUE
 123          SIGMA(KMAX) = SIGMA(J)
 124          WA(KMAX) = WA(J)
 125          K = IPVT(J)
 126          IPVT(J) = IPVT(KMAX)
 127          IPVT(KMAX) = K
 128    40    CONTINUE
 129 C
 130 C        COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE
 131 C        J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.
 132 C
 133          AJNORM = ENORM(M-J+1,A(J,J))
 134          IF (AJNORM .EQ. ZERO) GO TO 100
 135          IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM
 136          DO 50 I = J, M
 137             A(I,J) = A(I,J)/AJNORM
 138    50       CONTINUE
 139          A(J,J) = A(J,J) + ONE
 140 C
 141 C        APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
 142 C        AND UPDATE THE NORMS.
 143 C
 144          JP1 = J + 1
 145          IF (N .LT. JP1) GO TO 100
 146          DO 90 K = JP1, N
 147             SUM = ZERO
 148             DO 60 I = J, M
 149                SUM = SUM + A(I,J)*A(I,K)
 150    60          CONTINUE
 151             TEMP = SUM/A(J,J)
 152             DO 70 I = J, M
 153                A(I,K) = A(I,K) - TEMP*A(I,J)
 154    70          CONTINUE
 155             IF (.NOT.PIVOT .OR. SIGMA(K) .EQ. ZERO) GO TO 80
 156             TEMP = A(J,K)/SIGMA(K)
 157             SIGMA(K) = SIGMA(K)*SQRT(MAX(ZERO,ONE-TEMP**2))
 158             IF (P05*(SIGMA(K)/WA(K))**2 .GT. EPSMCH) GO TO 80
 159             SIGMA(K) = ENORM(M-J,A(JP1,K))
 160             WA(K) = SIGMA(K)
 161    80       CONTINUE
 162    90       CONTINUE
 163   100    CONTINUE
 164          SIGMA(J) = -AJNORM
 165   110    CONTINUE
 166       RETURN
 167 C
 168 C     LAST CARD OF SUBROUTINE QRFAC.
 169 C
 170       END