share/odepack/fortran/dheqr.f

   1 *DECK DHEQR
   2       SUBROUTINE DHEQR (A, LDA, N, Q, INFO, IJOB)
   3       INTEGER LDA, N, INFO, IJOB
   4       DOUBLE PRECISION A(LDA,*), Q(*)
   5 C-----------------------------------------------------------------------
   6 C     This routine performs a QR decomposition of an upper
   7 C     Hessenberg matrix A.  There are two options available:
   8 C
   9 C          (1)  performing a fresh decomposition
  10 C          (2)  updating the QR factors by adding a row and a
  11 C               column to the matrix A.
  12 C-----------------------------------------------------------------------
  13 C     DHEQR decomposes an upper Hessenberg matrix by using Givens
  14 C     rotations.
  15 C
  16 C     On entry
  17 C
  18 C        A       DOUBLE PRECISION(LDA, N)
  19 C                the matrix to be decomposed.
  20 C
  21 C        LDA     INTEGER
  22 C                the leading dimension of the array  A .
  23 C
  24 C        N       INTEGER
  25 C                A is an (N+1) by N Hessenberg matrix.
  26 C
  27 C        IJOB    INTEGER
  28 C                = 1     means that a fresh decomposition of the
  29 C                        matrix A is desired.
  30 C                .ge. 2  means that the current decomposition of A
  31 C                        will be updated by the addition of a row
  32 C                        and a column.
  33 C     On return
  34 C
  35 C        A       the upper triangular matrix R.
  36 C                The factorization can be written Q*A = R, where
  37 C                Q is a product of Givens rotations and R is upper
  38 C                triangular.
  39 C
  40 C        Q       DOUBLE PRECISION(2*N)
  41 C                the factors c and s of each Givens rotation used
  42 C                in decomposing A.
  43 C
  44 C        INFO    INTEGER
  45 C                = 0  normal value.
  46 C                = k  if  A(k,k) .eq. 0.0 .  This is not an error
  47 C                     condition for this subroutine, but it does
  48 C                     indicate that DHELS will divide by zero
  49 C                     if called.
  50 C
  51 C     Modification of LINPACK, by Peter Brown, LLNL.
  52 C     Written 1/13/86.  This version dated 6/20/01.
  53 C-----------------------------------------------------------------------
  54       INTEGER I, IQ, J, K, KM1, KP1, NM1
  55       DOUBLE PRECISION C, S, T, T1, T2
  56 C
  57       IF (IJOB .GT. 1) GO TO 70
  58 C
  59 C A new facorization is desired.
  60 C
  61 C     QR decomposition without pivoting
  62 C
  63       INFO = 0
  64       DO 60 K = 1, N
  65          KM1 = K - 1
  66          KP1 = K + 1
  67 C
  68 C           Compute kth column of R.
  69 C           First, multiply the kth column of A by the previous
  70 C           k-1 Givens rotations.
  71 C
  72             IF (KM1 .LT. 1) GO TO 20
  73             DO 10 J = 1, KM1
  74               I = 2*(J-1) + 1
  75               T1 = A(J,K)
  76               T2 = A(J+1,K)
  77               C = Q(I)
  78               S = Q(I+1)
  79               A(J,K) = C*T1 - S*T2
  80               A(J+1,K) = S*T1 + C*T2
  81    10         CONTINUE
  82 C
  83 C           Compute Givens components c and s
  84 C
  85    20       CONTINUE
  86             IQ = 2*KM1 + 1
  87             T1 = A(K,K)
  88             T2 = A(KP1,K)
  89             IF (T2 .NE. 0.0D0) GO TO 30
  90               C = 1.0D0
  91               S = 0.0D0
  92               GO TO 50
  93    30       CONTINUE
  94             IF (ABS(T2) .LT. ABS(T1)) GO TO 40
  95               T = T1/T2
  96               S = -1.0D0/SQRT(1.0D0+T*T)
  97               C = -S*T
  98               GO TO 50
  99    40       CONTINUE
 100               T = T2/T1
 101               C = 1.0D0/SQRT(1.0D0+T*T)
 102               S = -C*T
 103    50       CONTINUE
 104             Q(IQ) = C
 105             Q(IQ+1) = S
 106             A(K,K) = C*T1 - S*T2
 107             IF (A(K,K) .EQ. 0.0D0) INFO = K
 108    60 CONTINUE
 109       RETURN
 110 C
 111 C The old factorization of A will be updated.  A row and a column
 112 C has been added to the matrix A.
 113 C N by N-1 is now the old size of the matrix.
 114 C
 115   70  CONTINUE
 116       NM1 = N - 1
 117 C
 118 C Multiply the new column by the N previous Givens rotations.
 119 C
 120       DO 100 K = 1,NM1
 121         I = 2*(K-1) + 1
 122         T1 = A(K,N)
 123         T2 = A(K+1,N)
 124         C = Q(I)
 125         S = Q(I+1)
 126         A(K,N) = C*T1 - S*T2
 127         A(K+1,N) = S*T1 + C*T2
 128  100    CONTINUE
 129 C
 130 C Complete update of decomposition by forming last Givens rotation,
 131 C and multiplying it times the column vector (A(N,N), A(N+1,N)).
 132 C
 133       INFO = 0
 134       T1 = A(N,N)
 135       T2 = A(N+1,N)
 136       IF (T2 .NE. 0.0D0) GO TO 110
 137         C = 1.0D0
 138         S = 0.0D0
 139         GO TO 130
 140  110  CONTINUE
 141       IF (ABS(T2) .LT. ABS(T1)) GO TO 120
 142         T = T1/T2
 143         S = -1.0D0/SQRT(1.0D0+T*T)
 144         C = -S*T
 145         GO TO 130
 146  120  CONTINUE
 147         T = T2/T1
 148         C = 1.0D0/SQRT(1.0D0+T*T)
 149         S = -C*T
 150  130  CONTINUE
 151       IQ = 2*N - 1
 152       Q(IQ) = C
 153       Q(IQ+1) = S
 154       A(N,N) = C*T1 - S*T2
 155       IF (A(N,N) .EQ. 0.0D0) INFO = N
 156       RETURN
 157 C----------------------- End of Subroutine DHEQR -----------------------
 158       END