dheqr.f

Go to the documentation of this file.

C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number
C W-7405-Eng-48.
C
      SUBROUTINE dheqr (A, LDA, N, Q, INFO, IJOB)
C
C***BEGIN PROLOGUE  DHEQR
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     This routine performs a QR decomposition of an upper
C     Hessenberg matrix A.  There are two options available:
C
C          (1)  performing a fresh decomposition
C          (2)  updating the QR factors by adding a row and A
C               column to the matrix A.
C
C     DHEQR decomposes an upper Hessenberg matrix by using Givens
C     rotations.
C
C     On entry
C
C        A       DOUBLE PRECISION(LDA, N)
C                The matrix to be decomposed.
C
C        LDA     INTEGER
C                The leading dimension of the array A.
C
C        N       INTEGER
C                A is an (N+1) by N Hessenberg matrix.
C
C        IJOB    INTEGER
C                = 1     Means that a fresh decomposition of the
C                        matrix A is desired.
C                .GE. 2  Means that the current decomposition of A
C                        will be updated by the addition of a row
C                        and a column.
C     On return
C
C        A       The upper triangular matrix R.
C                The factorization can be written Q*A = R, where
C                Q is a product of Givens rotations and R is upper
C                triangular.
C
C        Q       DOUBLE PRECISION(2*N)
C                The factors C and S of each Givens rotation used
C                in decomposing A.
C
C        INFO    INTEGER
C                = 0  normal value.
C                = K  If  A(K,K) .EQ. 0.0.  This is not an error
C                     condition for this subroutine, but it does
C                     indicate that DHELS will divide by zero
C                     if called.
C
C     Modification of LINPACK.
C     Peter Brown, Lawrence Livermore Natl. Lab.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED (NONE)
C
C***END PROLOGUE  DHEQR
C
      INTEGER LDA, N, INFO, IJOB
      DOUBLE PRECISION A(LDA,*), Q(*)
      INTEGER I, IQ, J, K, KM1, KP1, NM1
      DOUBLE PRECISION C, S, T, T1, T2
C
      IF (ijob .GT. 1) GO TO 70
C-----------------------------------------------------------------------
C A new factorization is desired.
C-----------------------------------------------------------------------
C
C     QR decomposition without pivoting.
C
      info = 0
      DO 60 k = 1, n
         km1 = k - 1
         kp1 = k + 1
C
C           Compute Kth column of R.
C           First, multiply the Kth column of A by the previous
C           K-1 Givens rotations.
C
            IF (km1 .LT. 1) GO TO 20
            DO 10 j = 1, km1
              i = 2*(j-1) + 1
              t1 = a(j,k)
              t2 = a(j+1,k)
              c = q(i)
              s = q(i+1)
              a(j,k) = c*t1 - s*t2
              a(j+1,k) = s*t1 + c*t2
   10         CONTINUE
C
C           Compute Givens components C and S.
C
   20       CONTINUE
            iq = 2*km1 + 1
            t1 = a(k,k)
            t2 = a(kp1,k)
            IF (t2 .NE. 0.0d0) GO TO 30
              c = 1.0d0
              s = 0.0d0
              GO TO 50
   30       CONTINUE
            IF (abs(t2) .LT. abs(t1)) GO TO 40
              t = t1/t2
              s = -1.0d0/sqrt(1.0d0+t*t)
              c = -s*t
              GO TO 50
   40       CONTINUE
              t = t2/t1
              c = 1.0d0/sqrt(1.0d0+t*t)
              s = -c*t
   50       CONTINUE
            q(iq) = c
            q(iq+1) = s
            a(k,k) = c*t1 - s*t2
            IF (a(k,k) .EQ. 0.0d0) info = k
   60 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C The old factorization of A will be updated.  A row and a column
C has been added to the matrix A.
C N by N-1 is now the old size of the matrix.
C-----------------------------------------------------------------------
  70  CONTINUE
      nm1 = n - 1
C-----------------------------------------------------------------------
C Multiply the new column by the N previous Givens rotations.
C-----------------------------------------------------------------------
      DO 100 k = 1,nm1
        i = 2*(k-1) + 1
        t1 = a(k,n)
        t2 = a(k+1,n)
        c = q(i)
        s = q(i+1)
        a(k,n) = c*t1 - s*t2
        a(k+1,n) = s*t1 + c*t2
 100    CONTINUE
C-----------------------------------------------------------------------
C Complete update of decomposition by forming last Givens rotation,
C and multiplying it times the column vector (A(N,N),A(NP1,N)).
C-----------------------------------------------------------------------
      info = 0
      t1 = a(n,n)
      t2 = a(n+1,n)
      IF (t2 .NE. 0.0d0) GO TO 110
        c = 1.0d0
        s = 0.0d0
        GO TO 130
 110  CONTINUE
      IF (abs(t2) .LT. abs(t1)) GO TO 120
        t = t1/t2
        s = -1.0d0/sqrt(1.0d0+t*t)
        c = -s*t
        GO TO 130
 120  CONTINUE
        t = t2/t1
        c = 1.0d0/sqrt(1.0d0+t*t)
        s = -c*t
 130  CONTINUE
      iq = 2*n - 1
      q(iq) = c
      q(iq+1) = s
      a(n,n) = c*t1 - s*t2
      IF (a(n,n) .EQ. 0.0d0) info = n
      RETURN
C
C------END OF SUBROUTINE DHEQR------------------------------------------
      END