This routine performs a QR decomposition of an upper Hessenberg matrix a
using Givens
rotations. There are two options available:
a
.Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=rk), | intent(inout) | :: | a(lda,n) |
On entry, the Hessenberg matrix to be decomposed. On return, the upper triangular
matrix R from the QR decomposition of |
||
integer, | intent(in) | :: | lda |
Leading dimension of |
||
integer, | intent(in) | :: | n |
Number of columns in |
||
real(kind=rk), | intent(out) | :: | q(2*n) |
Coefficients of the Givens rotations used in decomposing |
||
integer, | intent(out) | :: | info |
Info flag.
|
||
integer, | intent(in) | :: | ijob |
Job flag.
|
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public | :: | i | ||||
integer, | public | :: | iq | ||||
integer, | public | :: | j | ||||
integer, | public | :: | k | ||||
integer, | public | :: | km1 | ||||
integer, | public | :: | kp1 | ||||
integer, | public | :: | nm1 | ||||
real(kind=rk), | public | :: | c | ||||
real(kind=rk), | public | :: | s | ||||
real(kind=rk), | public | :: | t | ||||
real(kind=rk), | public | :: | t1 | ||||
real(kind=rk), | public | :: | t2 |
pure subroutine dheqr(a, lda, n, q, info, ijob) !! This routine performs a QR decomposition of an upper Hessenberg matrix `a` using Givens !! rotations. 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`. use daskr_kinds, only: rk, zero, one implicit none real(rk), intent(inout) :: a(lda, n) !! On entry, the Hessenberg matrix to be decomposed. On return, the upper triangular !! matrix R from the QR decomposition of `a`. integer, intent(in) :: lda !! Leading dimension of `a`. integer, intent(in) :: n !! Number of columns in `a` (originally a matrix with shape `(n+1, n)`). real(rk), intent(out) :: q(2*n) !! Coefficients of the Givens rotations used in decomposing `a`. integer, intent(out) :: info !! Info flag. !! `info = 0`, if successful. !! `info = k`, if `a(k,k) = 0`; this is not an error condition for this subroutine, but !! it does indicate that [[dhels]] will divide by zero if called. integer, intent(in) :: ijob !! Job flag. !! `ijob = 1`, means that a fresh decomposition of the matrix `a` is desired. !! `ijob >= 2`, means that the current decomposition of `a` will be updated by the !! addition of a row and a column. integer :: i, iq, j, k, km1, kp1, nm1 real(rk) :: c, s, t, t1, t2 ! A new factorization is desired. if (ijob == 1) then ! 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) goto 20 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 end do ! Compute Givens components C and S. 20 continue iq = 2*km1 + 1 t1 = a(k, k) t2 = a(kp1, k) if (t2 /= zero) goto 30 c = one s = zero goto 50 30 continue if (abs(t2) < abs(t1)) goto 40 t = t1/t2 s = -one/sqrt(one + t*t) c = -s*t goto 50 40 continue t = t2/t1 c = one/sqrt(one + 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) == zero) info = k end do ! 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. else if (ijob >= 2) then 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 end do ! Complete update of decomposition by forming last Givens rotation, ! and multiplying it times the column vector (A(N,N),A(NP1,N)). info = 0 t1 = a(n, n) t2 = a(n + 1, n) if (t2 /= zero) goto 110 c = one s = zero goto 130 110 continue if (abs(t2) < abs(t1)) goto 120 t = t1/t2 s = -one/sqrt(one + t*t) c = -s*t goto 130 120 continue t = t2/t1 c = one/sqrt(one + 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) == zero) info = n end if end subroutine dheqr