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

DHGEQZ Example Program Data
   5                                                  :Value of N
1.00        1.00        1.00        1.00        1.00
2.00        4.00        8.00       16.00       32.00
3.00        9.00       27.00       81.00      243.00
4.00       16.00       64.00      256.00     1024.00
5.00       25.00      125.00      625.00     3125.00  :End of matrix A
1.00        2.00        3.00        4.00        5.00
1.00        4.00        9.00       16.00       25.00
1.00        8.00       27.00       64.00      125.00
1.00       16.00       81.00      256.00      625.00
1.00       32.00      243.00     1024.00     3125.00  :End of matrix B

 DHGEQZ Example Program Results
 Matrix A after balancing
             1          2          3          4          5
 1      1.0000     1.0000     0.1000     0.1000     0.1000
 2      2.0000     4.0000     0.8000     1.6000     3.2000
 3      0.3000     0.9000     0.2700     0.8100     2.4300
 4      0.4000     1.6000     0.6400     2.5600    10.2400
 5      0.5000     2.5000     1.2500     6.2500    31.2500

 Matrix B after balancing
             1          2          3          4          5
 1      1.0000     2.0000     0.3000     0.4000     0.5000
 2      1.0000     4.0000     0.9000     1.6000     2.5000
 3      0.1000     0.8000     0.2700     0.6400     1.2500
 4      0.1000     1.6000     0.8100     2.5600     6.2500
 5      0.1000     3.2000     2.4300    10.2400    31.2500

 Matrix A in Hessenberg form
             1          2          3          4          5
 1     -2.1898    -0.3181     2.0547     4.7371    -4.6249
 2     -0.8395    -0.0426     1.7132     7.5194   -17.1850
 3      0.0000    -0.2846    -1.0101    -7.5927    26.4499
 4      0.0000     0.0000     0.0376     1.4070    -3.3643
 5      0.0000     0.0000     0.0000     0.3813    -0.9937

 Matrix B is triangular
             1          2          3          4          5
 1     -1.4248    -0.3476     2.1175     5.5813    -3.9269
 2      0.0000    -0.0782     0.1189     8.0940   -15.2928
 3      0.0000     0.0000     1.0021   -10.9356    26.5971
 4      0.0000     0.0000     0.0000     0.5820    -0.0730
 5      0.0000     0.0000     0.0000     0.0000     0.5321

 Minimal required LWORK =      5
 Actual value of  LWORK =     30

 Generalized eigenvalues
    1     ( -2.437,  0.000)
    2     (  0.607,  0.795)
    3     (  0.607, -0.795)
    4     (  1.000,  0.000)
    5     ( -0.410,  0.000)

    Program dhgeqz_example

!     DHGEQZ Example Program Text

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

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen
      Use lapack_interfaces, Only: dgeqrf, dggbal, dgghd3, dhgeqz, dormqr
      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
      Character (1) :: compq, compz, job
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
        beta(:), lscale(:), q(:, :), rscale(:), tau(:), work(:), z(:, :)
!     .. Intrinsic Procedures ..
      Intrinsic :: nint
!     .. Executable Statements ..
      Write (nout, *) 'DHGEQZ 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 (alphai(n), alphar(n), beta(n), a(lda,n), lscale(n), &
        q(ldq,ldq), rscale(n), b(ldb,n), tau(n), work(lwork), z(ldz,ldz))

!     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 dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, &
        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_real_gen('General', ' ', n, n, a, lda, &
        'Matrix A after balancing', ifail)

      Write (nout, *)
      Flush (nout)

!     Matrix B after balancing

      ifail = 0
      Call nagf_file_print_matrix_real_gen('General', ' ', n, n, b, ldb, &
        'Matrix B after balancing', ifail)

      Write (nout, *)
      Flush (nout)

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

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

!     Apply the orthogonal transformation to matrix A
      Call dormqr('L', 'T', 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 dgghd3(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_real_gen('General', ' ', n, n, a, lda, &
        'Matrix A in Hessenberg form', ifail)

      Write (nout, *)
      Flush (nout)

!     Matrix B (T) in generalized Hessenberg form.

      ifail = 0
      Call nagf_file_print_matrix_real_gen('General', ' ', n, n, b, ldb, &
        'Matrix B is triangular', ifail)

!     Routine DHGEQZ
!     Workspace query: jwork = -1

      jwork = -1
      job = 'E'

      Call dhgeqz(job, compq, compz, n, ilo, ihi, a, lda, b, ldb, alphar, &
        alphai, beta, q, ldq, z, ldz, work, jwork, info)

      Write (nout, *)
      Write (nout, 100) nint(work(1))
      Write (nout, 110) lwork
      Write (nout, *)

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

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

        Call dhgeqz(job, compq, compz, n, ilo, ihi, a, lda, b, ldb, alphar, &
          alphai, beta, q, ldq, z, ldz, work, lwork, info)

!       Print the generalized eigenvalues

        Write (nout, 120)

        Do i = 1, n
          If (beta(i)/=0.0E0_dp) Then
            Write (nout, 130) i, '(', alphar(i)/beta(i), ',', &
              alphai(i)/beta(i), ')'
          Else
            Write (nout, 150) i
          End If
        End Do
      Else
        Write (nout, 140)
      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, A, F7.3, A, F7.3, A)
140   Format (1X, 'Insufficient workspace allocated for call to DHGEQZ')
150   Format (1X, I4, 'Eigenvalue is infinite')
    End Program