share/hompack/fortran/qrfaqf.f

   1         SUBROUTINE QRFAQF(QT,R,N,IFLAG)
   2 C
   3 C SUBROUTINE  QRFAQF  COMPUTES THE QR FACTORIZATION OF A MATRIX  A,
   4 C WHERE  R  IS AN UPPER TRIANGULAR MATRIX, AND  Q  IS AN ORTHOGONAL
   5 C MATRIX WHICH IS THE PRODUCT OF N-1 HOUSEHOLDER TRANSFORMATIONS
   6 C
   7 C   Q=H1*H2*...*H(N-1).
   8 C
   9 C THE ROUTINE HAS TWO MAJOR STEPS.  FIRST, THE QR FACTORIZATION
  10 C OF  A  IS COMPUTED, RESULTING IN DEFINING THE VECTOR  R,  AND
  11 C STORING INFORMATION IN THE LOWER TRIANGLE OF  QT  WHICH WILL
  12 C ENABLE THE CONSTRUCTION OF  Q TRANSPOSE.
  13 C
  14 C THE SECOND STEP CONSTRUCTS  Q TRANSPOSE FROM THE INFORMATION
  15 C STORED IN  QT,  AND PLACES IT IN  QT.
  16 C
  17 C THE INFORMATION STORED IN THE LOWER TRIANGLE OF  QT  DURING THE FIRST
  18 C STEP ARE THE VECTORS  UJ, WHICH DEFINE THE HOUSEHOLDER TRANSFORMATIONS
  19 C
  20 C                   T
  21 C   HJ = I - (UJ*UJ  / PJ),  WHERE UJ[I]=0 FOR I=1...J-1,
  22 C                             UJ[I]=QT[I,J], FOR I=J...N,
  23 C                             PJ = THE JTH COMPONENT OF UJ.
  24 C
  25 C
  26 C ON INPUT:
  27 C
  28 C QT(1:N,1:N)  CONTAINS THE MATRIX  A  TO BE FACTORED.
  29 C
  30 C R(1:N*(N+1)/2)  IS UNDEFINED.
  31 C
  32 C N  IS THE DIMENSION OF THE MATRIX TO BE FACTORED.
  33 C
  34 C IFLAG  IS UNDEFINED.
  35 C
  36 C
  37 C ON OUTPUT:
  38 C
  39 C QT  CONTAINS  Q TRANSPOSE.
  40 C
  41 C R(1:N*(N+1)/2)  CONTAINS THE UPPER TRIANGLE OF  R  STORED BY ROWS.
  42 C
  43 C N  IS UNCHANGED.
  44 C
  45 C IFLAG = 4  IF THE MATRIX A IS SINGULAR.  OTHERWISE,  IFLAG
  46 C    IS UNCHANGED.
  47 C
  48 C
  49 C CALLS  DAXPY, DCOPY, DDOT, DNRM2, DSCAL.
  50 C
  51 C ***** DECLARATIONS *****
  52 C
  53 C     FUNCTION DECLARATIONS
  54 C
  55         DOUBLE PRECISION DDOT, DNRM2
  56 C
  57 C     LOCAL VARIABLES
  58 C
  59         DOUBLE PRECISION ONE, TAU, TEMP
  60         INTEGER I, J, K, INDEXR, ISIGN
  61 C
  62 C     SCALAR ARGUMENTS
  63 C
  64         INTEGER N, IFLAG
  65 C
  66 C     ARRAY DECLARATIONS
  67 C
  68         DOUBLE PRECISION QT(N,N),R(N)
  69 C
  70 C ***** END OF DECLARATIONS *****
  71 C
  72 C ***** FIRST EXECUTABLE STATEMENT *****
  73 C
  74         ONE = 1.0
  75 C
  76 C ***** CALCULATION OF QR DECOMPOSITION, PLACING R IN THE VECTOR *****
  77 C       R, AND PLACING THE UJ VECTORS IN THE LOWER TRIANGLE OF
  78 C       QT.
  79 C
  80         INDEXR = 1
  81         DO 20 K=1,N-1
  82           TEMP = DNRM2(N-K+1,QT(K,K),1)
  83           IF (TEMP .EQ. 0.0) THEN
  84 C
  85 C           MATRIX IS SINGULAR, SET  IFLAG  AND RETURN.
  86 C
  87             IFLAG = 4
  88             RETURN
  89           ELSE
  90 C
  91 C           FORM QK AND PREMULTIPLY  QT  BY IT.
  92 C                                                         T
  93 C           UK = EK - ISIGN*X/||X||,  WHERE  HK = I-(UK*UK /PK),
  94 C              PK = THE KTH COMPONENT OF UK,
  95 C              EK = THE KTH NATURAL BASIS VECTOR,
  96 C              X  = THE KTH COLUMN OF THE MATRIX  H(K-1)...H2*H1*QT,
  97 C              ISIGN = THE SIGN OF PK.
  98 C
  99 C           GET SIGN.
 100 C
 101             ISIGN = SIGN(ONE,QT(K,K))
 102 C
 103 C           COMPUTE R(K,K).
 104 C
 105             R(INDEXR) = -ISIGN*TEMP
 106 C
 107 C           UPDATE KTH COLUMN.
 108 C
 109             TEMP = ISIGN/TEMP
 110             CALL DSCAL(N-K+1,TEMP,QT(K,K),1)
 111             QT(K,K) = QT(K,K) + 1.0
 112 C
 113 C           UPDATE THE K+1ST - NTH COLUMNS OF  QT, AND  R.
 114 C
 115             INDEXR = INDEXR + 1
 116             DO 10 J=K+1,N
 117               TAU = DDOT(N-K+1,QT(K,K),1,QT(K,J),1)/QT(K,K)
 118               R(INDEXR) = QT(K,J) - TAU*QT(K,K)
 119               INDEXR = INDEXR + 1
 120               CALL DAXPY(N-K,-TAU,QT(K+1,K),1,QT(K+1,J),1)
 121   10        CONTINUE
 122           END IF
 123   20    CONTINUE
 124         IF (QT(N,N) .EQ. 0.0) THEN
 125 C
 126 C         MATRIX IS SINGULAR, SET  IFLAG  AND RETURN.
 127 C
 128           IFLAG = 4
 129           RETURN
 130         END IF
 131         R(INDEXR) = QT(N,N)
 132 C
 133 C ***** END OF FACTORING STEP *****
 134 C
 135 C ***** CONSTRUCT Q TRANSPOSE IN QT *****
 136 C
 137 C FORM Q BY MULTIPLYING ((I*H(N-1))*...)*H1.
 138 C THIS IS DONE IN PLACE IN  QT  BY UPDATING ONLY THE LOWER
 139 C RIGHT HAND CORNER OF QT  (QT(K,K) TO QT(N,N)).
 140 C
 141 C
 142         QT(N,N) = 1.0
 143         DO 40 K=N-1,1,-1
 144 C
 145 C         MULTIPLY  QT  BY H(K).
 146 C
 147           TEMP = QT(K,K)
 148 C
 149 C         UPDATE ROW K.
 150 C
 151           QT(K,K) = 1.0-QT(K,K)
 152           CALL DCOPY(N-K,QT(K+1,K),1,QT(K,K+1),N)
 153           CALL DSCAL(N-K,-ONE,QT(K,K+1),N)
 154 C
 155 C         UPDATE REMAINING ROWS.
 156 C
 157           DO 30 I=N,K+1,-1
 158             TAU = -DDOT(N-K,QT(I,K+1),N,QT(K,K+1),N)
 159             QT(I,K) = -TAU
 160             TAU = TAU/TEMP
 161             CALL DAXPY(N-K,TAU,QT(K,K+1),N,QT(I,K+1),N)
 162   30      CONTINUE
 163   40    CONTINUE
 164 C
 165 C ***** END OF Q TRANSPOSE CONSTRUCTION *****
 166 C
 167         RETURN
 168 C
 169 C ***** END OF SUBROUTINE QRFAQF *****
 170         END