計算ルーチン: ２つの複素行列から縮約された一般上複素ヘッセンベルグの固有値とシュール分解

ZHGEQZ Example Program Data
   4                                                    : n

(  1.0, 3.0)  (  1.0, 4.0)  (  1.0, 5.0)  (  1.0, 6.0)
(  2.0, 2.0)  (  4.0, 3.0)  (  8.0, 4.0)  ( 16.0, 5.0)
(  3.0, 1.0)  (  9.0, 2.0)  ( 27.0, 3.0)  ( 81.0, 4.0)
(  4.0, 0.0)  ( 16.0, 1.0)  ( 64.0, 2.0)  (256.0, 3.0)  : A

(  1.0, 0.0)  (  2.0, 1.0)  (  3.0, 2.0)  (  4.0, 3.0)
(  1.0, 1.0)  (  4.0, 2.0)  (  9.0, 3.0)  ( 16.0, 4.0)
(  1.0, 2.0)  (  8.0, 3.0)  ( 27.0, 4.0)  ( 64.0, 5.0)
(  1.0, 3.0)  ( 16.0, 4.0)  ( 81.0, 5.0)  (256.0, 6.0)  : B

 ZHGEQZ Example Program Results
 Matrix A after balancing
                    1                 2                 3                 4
 1  ( 1.0000, 3.0000) ( 1.0000, 4.0000) ( 0.1000, 0.5000) ( 0.1000, 0.6000)
 2  ( 2.0000, 2.0000) ( 4.0000, 3.0000) ( 0.8000, 0.4000) ( 1.6000, 0.5000)
 3  ( 0.3000, 0.1000) ( 0.9000, 0.2000) ( 0.2700, 0.0300) ( 0.8100, 0.0400)
 4  ( 0.4000, 0.0000) ( 1.6000, 0.1000) ( 0.6400, 0.0200) ( 2.5600, 0.0300)

 Matrix B after balancing
                    1                 2                 3                 4
 1  ( 1.0000, 0.0000) ( 2.0000, 1.0000) ( 0.3000, 0.2000) ( 0.4000, 0.3000)
 2  ( 1.0000, 1.0000) ( 4.0000, 2.0000) ( 0.9000, 0.3000) ( 1.6000, 0.4000)
 3  ( 0.1000, 0.2000) ( 0.8000, 0.3000) ( 0.2700, 0.0400) ( 0.6400, 0.0500)
 4  ( 0.1000, 0.3000) ( 1.6000, 0.4000) ( 0.8100, 0.0500) ( 2.5600, 0.0600)

 Matrix A in Hessenberg form
                    1                 2                 3                 4
 1  ( -2.868, -1.595) ( -0.809, -0.328) ( -4.900, -0.987) ( -0.048,  1.163)
 2  ( -2.672,  0.595) ( -0.790,  0.049) ( -4.955, -0.163) ( -0.439, -0.574)
 3  (  0.000,  0.000) ( -0.098, -0.011) ( -1.168, -0.137) ( -1.756, -0.205)
 4  (  0.000,  0.000) (  0.000,  0.000) (  0.087,  0.004) (  0.032,  0.001)

 Matrix B is triangular
                    1                 2                 3                 4
 1  ( -1.775,  0.000) ( -0.721,  0.043) ( -5.021,  1.190) ( -0.145,  0.726)
 2  (  0.000,  0.000) ( -0.218,  0.035) ( -2.541, -0.146) ( -0.823, -0.418)
 3  (  0.000,  0.000) (  0.000,  0.000) ( -1.396, -0.163) ( -1.747, -0.204)
 4  (  0.000,  0.000) (  0.000,  0.000) (  0.000,  0.000) ( -0.100, -0.004)

 Minimal required LWORK =      4
 Actual value of  LWORK =     24

 Generalized eigenvalues

    1     ( -0.635,  1.653)
    2     (  0.493,  0.910)
    3     (  0.458, -0.843)
    4     (  0.674, -0.050)

    Program zhgeqz_example

!     ZHGEQZ Example Program Text

!     Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_sort_realvec_rank, &
        nagf_file_print_matrix_complex_gen_comp, &
        nagf_sort_cmplxvec_rank_rearrange
      Use lapack_interfaces, Only: zgeqrf, zggbal, zgghd3, zhgeqz, zunmqr
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, ifail, ihi, ilo, info, irows, jwork, lda, ldb, ldq, ldz, &
        lwork, n, ni
      Character (1) :: compq, compz, job
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), alpha(:), b(:, :), beta(:), &
        e(:), q(:, :), tau(:), work(:), z(:, :)
      Real (Kind=dp), Allocatable :: emod(:), lscale(:), rscale(:), rwork(:)
      Integer, Allocatable :: irank(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, aimag, nint, real
!     .. Executable Statements ..
      Write (nout, *) 'ZHGEQZ Example Program Results'
      Flush (nout)

!     Skip heading in data file

      Read (nin, *)
      Read (nin, *) n
      ldq = 1
      ldz = 1
      lda = n
      ldb = n
      lwork = 6*n
      Allocate (a(lda,n), alpha(n), b(ldb,n), beta(n), q(ldq,ldq), tau(n), &
        work(lwork), z(ldz,ldz), lscale(n), rscale(n), rwork(6*n))

!     READ matrix A from data file
      Read (nin, *)(a(i,1:n), i=1, n)

!     READ matrix B from data file
      Read (nin, *)(b(i,1:n), i=1, n)

!     Balance matrix pair (A,B)
      job = 'B'
      Call zggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, rwork, &
        info)

!     Matrix A after balancing

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, n, a, &
        lda, 'Bracketed', 'F7.4', 'Matrix A after balancing', 'Integer', &
        rlabs, 'Integer', clabs, 80, 0, ifail)

      Write (nout, *)
      Flush (nout)

!     Matrix B after balancing

      ifail = 0
      Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, n, b, &
        ldb, 'Bracketed', 'F7.4', 'Matrix B after balancing', 'Integer', &
        rlabs, 'Integer', clabs, 80, 0, ifail)

      Write (nout, *)
      Flush (nout)

!     Reduce B to triangular form using QR
      irows = ihi + 1 - ilo

      Call zgeqrf(irows, irows, b(ilo,ilo), ldb, tau, work, lwork, info)

!     Apply the orthogonal transformation to A
      Call zunmqr('L', 'C', irows, irows, irows, b(ilo,ilo), ldb, tau, &
        a(ilo,ilo), lda, work, lwork, info)

!     Compute the generalized Hessenberg form of (A,B) -> (H,T)
      compq = 'N'
      compz = 'N'

      Call zgghd3(compq, compz, irows, 1, irows, a(ilo,ilo), lda, b(ilo,ilo), &
        ldb, q, ldq, z, ldz, work, lwork, info)

!     Matrix A (H) in generalized Hessenberg form
      ifail = 0
      Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, n, a, &
        lda, 'Bracketed', 'F7.3', 'Matrix A in Hessenberg form', 'Integer', &
        rlabs, 'Integer', clabs, 80, 0, ifail)

      Write (nout, *)
      Flush (nout)

!     Matrix B (T) in generalized Hessenberg form
      ifail = 0
      Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, n, b, &
        ldb, 'Bracketed', 'F7.3', 'Matrix B is triangular', 'Integer', rlabs, &
        'Integer', clabs, 80, 0, ifail)

!     Routine ZHGEQZ
!     Workspace query: jwork = -1

      jwork = -1
      job = 'E'
      Call zhgeqz(job, compq, compz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, &
        q, ldq, z, ldz, work, jwork, rwork, info)
      Write (nout, *)
      Write (nout, 100) nint(real(work(1)))
      Write (nout, 110) lwork
      Write (nout, *)
      Write (nout, 120)
      Write (nout, *)
      Flush (nout)

!     Compute the generalized Schur form
!     if the workspace lwork is adequate

      If (nint(real(work(1)))<=lwork) Then

        Call zhgeqz(job, compq, compz, n, ilo, ihi, a, lda, b, ldb, alpha, &
          beta, q, ldq, z, ldz, work, lwork, rwork, info)

!       Print the generalized eigenvalues in descending size order
!       Note: the actual values of beta are real and non-negative

!       Calculate the moduli of the finite eigenvalues.
        Allocate (e(n), emod(n), irank(n))
        ni = 0
        Do i = 1, n
          If (real(beta(i))/=0.0_dp) Then
            ni = ni + 1
            e(ni) = alpha(i)/beta(i)
            emod(ni) = abs(e(ni))
          Else
            Write (nout, 130) i
          End If
        End Do

!       Rearrange the finite eigenvalues in descending order of modulus.
        ifail = 0
        Call nagf_sort_realvec_rank(emod, 1, ni, 'Descending', irank, ifail)
        ifail = 0
        Call nagf_sort_cmplxvec_rank_rearrange(e, 1, ni, irank, ifail)

        Write (nout, 140)(i, '(', real(e(i)), ',', aimag(e(i)), ')', i=1, ni)
      Else
        Write (nout, 150)
      End If

100   Format (1X, 'Minimal required LWORK = ', I6)
110   Format (1X, 'Actual value of  LWORK = ', I6)
120   Format (1X, 'Generalized eigenvalues')
130   Format (1X, I4, 5X, 'Infinite eigenvalue')
140   Format (1X, I4, 5X, A, F7.3, A, F7.3, A)
150   Format (1X, 'Insufficient workspace allocated for call to ZHGEQZ')
    End Program