dheqr Subroutine

subroutine dheqr(A, Lda, N, Q, Info, Ijob)

This routine performs a QR decomposition of an upper Hessenberg matrix A. There are two options available:

      (1)  performing a fresh decomposition
      (2)  updating the QR factors by adding a row and a
           column to the matrix A.

DHEQR decomposes an upper Hessenberg matrix by using Givens rotations.

On entry

    A       DOUBLE PRECISION(LDA, N)
            the matrix to be decomposed.

    LDA     INTEGER
            the leading dimension of the array  A .

    N       INTEGER
            A is an (N+1) by N Hessenberg matrix.

    IJOB    INTEGER
            = 1     means that a fresh decomposition of the
                    matrix A is desired.
            .ge. 2  means that the current decomposition of A
                    will be updated by the addition of a row
                    and a column.

On return

    A       the upper triangular matrix R.
            The factorization can be written Q*A = R, where
            Q is a product of Givens rotations and R is upper
            triangular.

    Q       DOUBLE PRECISION(2*N)
            the factors c and s of each Givens rotation used
            in decomposing A.

    INFO    INTEGER
            = 0  normal value.
            = k  if  A(k,k) .eq. 0.0 .  This is not an error
                 condition for this subroutine, but it does
                 indicate that DHELS will divide by zero
                 if called.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(inout) :: A(Lda,*)
integer, intent(in) :: Lda
integer, intent(in) :: N
real(kind=dp), intent(inout) :: Q(*)
integer, intent(out) :: Info
integer, intent(in) :: Ijob

Variables

Type Visibility Attributes Name Initial
real(kind=dp), public :: c
integer, public :: i
integer, public :: iq
integer, public :: j
integer, public :: k
integer, public :: km1
integer, public :: kp1
integer, public :: nm1
real(kind=dp), public :: s
real(kind=dp), public :: t
real(kind=dp), public :: t1
real(kind=dp), public :: t2

Source Code

subroutine dheqr(A,Lda,N,Q,Info,Ijob)
!
integer,intent(in)          :: Lda
real(kind=dp),intent(inout) :: A(Lda,*)
integer,intent(in)          :: N
real(kind=dp),intent(inout) :: Q(*)
integer,intent(out)         :: Info
integer,intent(in)          :: Ijob
!
real(kind=dp) :: c , s , t , t1 , t2
integer :: i , iq , j , k , km1 , kp1 , nm1
!
if ( Ijob>1 ) then
!
!  The old factorization of A will be updated.  A row and a column
!  has been added to the matrix A.
!  N by N-1 is now the old size of the matrix.
!
   nm1 = N - 1
!
!  Multiply the new column by the N previous Givens rotations.
!
   do 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
   enddo
!
!  Complete update of decomposition by forming last Givens rotation,
!  and multiplying it times the column vector (A(N,N), A(N+1,N)).
!
   Info = 0
   t1 = A(N,N)
   t2 = A(N+1,N)
   if ( t2==0.0D0 ) then
      c = 1.0D0
      s = 0.0D0
   elseif ( abs(t2)<abs(t1) ) then
      t = t2/t1
      c = 1.0D0/sqrt(1.0D0+t*t)
      s = -c*t
   else
      t = t1/t2
      s = -1.0D0/sqrt(1.0D0+t*t)
      c = -s*t
   endif
   iq = 2*N - 1
   Q(iq) = c
   Q(iq+1) = s
   A(N,N) = c*t1 - s*t2
   if ( A(N,N)==0.0D0 ) Info = N
else
!
!  A new facorization is desired.
!
!      QR decomposition without pivoting
!
   Info = 0
   do k = 1 , N
      km1 = k - 1
      kp1 = k + 1
!
!            Compute kth column of R.
!            First, multiply the kth column of A by the previous
!            k-1 Givens rotations.
!
      if ( km1>=1 ) then
         do 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
         enddo
      endif
!
!            Compute Givens components c and s
!
      iq = 2*km1 + 1
      t1 = A(k,k)
      t2 = A(kp1,k)
      if ( t2==0.0D0 ) then
         c = 1.0D0
         s = 0.0D0
      elseif ( abs(t2)<abs(t1) ) then
         t = t2/t1
         c = 1.0D0/sqrt(1.0D0+t*t)
         s = -c*t
      else
         t = t1/t2
         s = -1.0D0/sqrt(1.0D0+t*t)
         c = -s*t
      endif
      Q(iq) = c
      Q(iq+1) = s
      A(k,k) = c*t1 - s*t2
      if ( A(k,k)==0.0D0 ) Info = k
   enddo
   return
endif
end subroutine dheqr