SUBROUTINE QRFAQF(QT,R,N,IFLAG)
C
C SUBROUTINE  QRFAQF  COMPUTES THE QR FACTORIZATION OF A MATRIX  A,
C WHERE  R  IS AN UPPER TRIANGULAR MATRIX, AND  Q  IS AN ORTHOGONAL
C MATRIX WHICH IS THE PRODUCT OF N-1 HOUSEHOLDER TRANSFORMATIONS
C
C   Q=H1*H2*...*H(N-1).
C
C THE ROUTINE HAS TWO MAJOR STEPS.  FIRST, THE QR FACTORIZATION
C OF  A  IS COMPUTED, RESULTING IN DEFINING THE VECTOR  R,  AND
C STORING INFORMATION IN THE LOWER TRIANGLE OF  QT  WHICH WILL
C ENABLE THE CONSTRUCTION OF  Q TRANSPOSE.
C
C THE SECOND STEP CONSTRUCTS  Q TRANSPOSE FROM THE INFORMATION
C STORED IN  QT,  AND PLACES IT IN  QT.
C
C THE INFORMATION STORED IN THE LOWER TRIANGLE OF  QT  DURING THE FIRST
C STEP ARE THE VECTORS  UJ, WHICH DEFINE THE HOUSEHOLDER TRANSFORMATIONS
C
C                   T
C   HJ = I - (UJ*UJ  / PJ),  WHERE UJ[I]=0 FOR I=1...J-1,
C                             UJ[I]=QT[I,J], FOR I=J...N,
C                             PJ = THE JTH COMPONENT OF UJ.
C
C
C ON INPUT:
C
C QT(1:N,1:N)  CONTAINS THE MATRIX  A  TO BE FACTORED.
C
C R(1:N*(N+1)/2)  IS UNDEFINED.
C
C N  IS THE DIMENSION OF THE MATRIX TO BE FACTORED.
C
C IFLAG  IS UNDEFINED.
C
C
C ON OUTPUT:
C
C QT  CONTAINS  Q TRANSPOSE.
C
C R(1:N*(N+1)/2)  CONTAINS THE UPPER TRIANGLE OF  R  STORED BY ROWS.
C
C N  IS UNCHANGED.
C
C IFLAG = 4  IF THE MATRIX A IS SINGULAR.  OTHERWISE,  IFLAG
C    IS UNCHANGED.
C
C
C CALLS  DAXPY, DCOPY, DDOT, DNRM2, DSCAL.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION DDOT, DNRM2
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION ONE, TAU, TEMP
        INTEGER I, J, K, INDEXR, ISIGN
C
C     SCALAR ARGUMENTS
C
        INTEGER N, IFLAG
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION QT(N,N),R(N)
C
C ***** END OF DECLARATIONS *****
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
        ONE = 1.0
C
C ***** CALCULATION OF QR DECOMPOSITION, PLACING R IN THE VECTOR *****
C       R, AND PLACING THE UJ VECTORS IN THE LOWER TRIANGLE OF
C       QT.
C
        INDEXR = 1
        DO 20 K=1,N-1
          TEMP = DNRM2(N-K+1,QT(K,K),1)
          IF (TEMP .EQ. 0.0) THEN
C
C           MATRIX IS SINGULAR, SET  IFLAG  AND RETURN.
C
            IFLAG = 4
            RETURN
          ELSE
C
C           FORM QK AND PREMULTIPLY  QT  BY IT.
C                                                         T
C           UK = EK - ISIGN*X/||X||,  WHERE  HK = I-(UK*UK /PK),
C              PK = THE KTH COMPONENT OF UK,
C              EK = THE KTH NATURAL BASIS VECTOR,
C              X  = THE KTH COLUMN OF THE MATRIX  H(K-1)...H2*H1*QT,
C              ISIGN = THE SIGN OF PK.
C
C           GET SIGN.
C
            ISIGN = SIGN(ONE,QT(K,K))
C
C           COMPUTE R(K,K).
C
            R(INDEXR) = -ISIGN*TEMP
C
C           UPDATE KTH COLUMN.
C
            TEMP = ISIGN/TEMP
            CALL DSCAL(N-K+1,TEMP,QT(K,K),1)
            QT(K,K) = QT(K,K) + 1.0
C
C           UPDATE THE K+1ST - NTH COLUMNS OF  QT, AND  R.
C
            INDEXR = INDEXR + 1
            DO 10 J=K+1,N
              TAU = DDOT(N-K+1,QT(K,K),1,QT(K,J),1)/QT(K,K)
              R(INDEXR) = QT(K,J) - TAU*QT(K,K)
              INDEXR = INDEXR + 1
              CALL DAXPY(N-K,-TAU,QT(K+1,K),1,QT(K+1,J),1)
  10        CONTINUE
          END IF
  20    CONTINUE
        IF (QT(N,N) .EQ. 0.0) THEN
C
C         MATRIX IS SINGULAR, SET  IFLAG  AND RETURN.
C
          IFLAG = 4
          RETURN
        END IF
        R(INDEXR) = QT(N,N)
C
C ***** END OF FACTORING STEP *****
C
C ***** CONSTRUCT Q TRANSPOSE IN QT *****
C
C FORM Q BY MULTIPLYING ((I*H(N-1))*...)*H1.
C THIS IS DONE IN PLACE IN  QT  BY UPDATING ONLY THE LOWER
C RIGHT HAND CORNER OF QT  (QT(K,K) TO QT(N,N)).
C
C
        QT(N,N) = 1.0
        DO 40 K=N-1,1,-1
C
C         MULTIPLY  QT  BY H(K).
C
          TEMP = QT(K,K)
C
C         UPDATE ROW K.
C
          QT(K,K) = 1.0-QT(K,K)
          CALL DCOPY(N-K,QT(K+1,K),1,QT(K,K+1),N)
          CALL DSCAL(N-K,-ONE,QT(K,K+1),N)
C
C         UPDATE REMAINING ROWS.
C
          DO 30 I=N,K+1,-1
            TAU = -DDOT(N-K,QT(I,K+1),N,QT(K,K+1),N)
            QT(I,K) = -TAU
            TAU = TAU/TEMP
            CALL DAXPY(N-K,TAU,QT(K,K+1),N,QT(I,K+1),N)
  30      CONTINUE
  40    CONTINUE
C
C ***** END OF Q TRANSPOSE CONSTRUCTION *****
C
        RETURN
C
C ***** END OF SUBROUTINE QRFAQF *****
        END

.