実一般化非対称固有値問題: 実非対称行列ペアの一般化固有値 : (オプションで左／右一般化固有ベクトル)

DGGEV3 Example Program Data
   4                     :Value of N
   3.9  12.5 -34.5  -0.5
   4.3  21.5 -47.5   7.5
   4.3  21.5 -43.5   3.5
   4.4  26.0 -46.0   6.0 :End of matrix A
   1.0   2.0  -3.0   1.0
   1.0   3.0  -5.0   4.0
   1.0   3.0  -4.0   3.0
   1.0   3.0  -4.0   4.0 :End of matrix B

 DGGEV3 Example Program Results
 Eigenvalues:
             1          2          3          4
 1      3.0000     3.0000     4.0000     2.0000
        4.0000    -4.0000     0.0000     0.0000

Warning: Floating underflow occurred
 Right Eigenvectors (columns):
          1       2       3       4
 1   0.7122  0.7122  1.0000  1.0000
     0.0000  0.0000  0.0000  0.0000
 
 2   0.1424  0.1424  0.0111  0.0057
     0.0000 -0.0000  0.0000  0.0000
 
 3   0.0855  0.0855 -0.0333  0.0629
    -0.1140  0.1140 -0.0000  0.0000
 
 4   0.0855  0.0855  0.1556  0.0629
    -0.1140  0.1140  0.0000  0.0000

    Program dggev3_example

!     DGGEV3 Example Program Text

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

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_sort_cmplxvec_rank_rearrange, &
        nagf_sort_realmat_rank_rows, nagf_file_print_matrix_complex_gen, &
        nagf_file_print_matrix_real_gen
      Use lapack_interfaces, Only: dggev3
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=dp), Parameter :: zero = 0.0_dp
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Complex (Kind=dp) :: scal
      Integer :: i, ifail, info, j, k, lda, ldb, ldvr, lwork, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: eigval(:), eigvec(:, :)
      Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
        beta(:), vr(:, :), work(:)
      Real (Kind=dp) :: dummy(1, 1)
      Integer, Allocatable :: irank(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, all, cmplx, conjg, epsilon, maxloc, nint
!     .. Executable Statements ..
      Write (nout, *) 'DGGEV3 Example Program Results'
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      ldb = n
      ldvr = n
      Allocate (a(lda,n), alphai(n), alphar(n), b(ldb,n), beta(n), vr(ldvr,n), &
        eigvec(n,n), eigval(n), irank(n))

!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call dggev3('No left vectors', 'Vectors (right)', n, a, lda, b, ldb, &
        alphar, alphai, beta, dummy, 1, vr, ldvr, dummy, lwork, info)

      lwork = nint(dummy(1,1))
      Allocate (work(lwork))

!     Read in the matrices A and B

      Read (nin, *)(a(i,1:n), i=1, n)
      Read (nin, *)(b(i,1:n), i=1, n)

!     Solve the generalized eigenvalue problem

      Call dggev3('No left vectors', 'Vectors (right)', n, a, lda, b, ldb, &
        alphar, alphai, beta, dummy, 1, vr, ldvr, work, lwork, info)
      If (info>0) Then
        Write (nout, *)
        Write (nout, 110) 'Failure in DGGEV3. INFO =', info
        Go To 100
      End If

!     Re-normalize the eigenvectors, largest absolute element real
      j = 0
      Do i = 1, n
        If (alphai(i)==zero) Then
          eigvec(1:n, i) = cmplx(vr(1:n,i), zero, kind=dp)
        Else If (j==0) Then
          eigvec(1:n, i) = cmplx(vr(1:n,i), vr(1:n,i+1), kind=dp)
          j = 1
        Else
          eigvec(1:n, i) = cmplx(vr(1:n,i-1), -vr(1:n,i), kind=dp)
          j = 0
        End If
        work(1:n) = abs(eigvec(1:n,i))
        k = maxloc(work(1:n), 1)
        scal = conjg(eigvec(k,i))/abs(eigvec(k,i))
        eigvec(1:n, i) = eigvec(1:n, i)*scal
      End Do

!     If eigenvalues are finite, order by descending absolute values
      If (all(abs(beta(1:n))>epsilon(1.0E0_dp))) Then
!       add small amount to alphai to distinguish conjugates
        alphai(1:n) = alphai(1:n) + epsilon(1.0E0_dp)*10.0_dp
        eigval(1:n) = cmplx(alphar(1:n), alphai(1:n), kind=dp)
        eigval(1:n) = eigval(1:n)/beta(1:n)
        work(1:n) = abs(eigval(1:n))
        ifail = 0
        Call nagf_sort_realmat_rank_rows(work, n, 1, n, 1, 1, 'Descending', &
          irank, ifail)
        Call nagf_sort_cmplxvec_rank_rearrange(eigval, 1, n, irank, ifail)

!       Print ordered eigenvalues
        ifail = 0
        Call nagf_file_print_matrix_complex_gen('Gen', ' ', 1, n, eigval, 1, &
          'Eigenvalues:', ifail)

!       Order the eigenvectors in the same way and print
        Do j = 1, n
          eigval(1:n) = eigvec(j, 1:n)
          Call nagf_sort_cmplxvec_rank_rearrange(eigval, 1, n, irank, ifail)
          eigvec(j, 1:n) = eigval(1:n)
        End Do

        Write (nout, *)
        Flush (nout)
        ifail = 0
        Call nagf_file_print_matrix_complex_gen('Gen', ' ', n, n, eigvec, n, &
          'Right Eigenvectors (columns):', ifail)
      Else
        Write (nout, *) 'Some of the eigenvalues are infinite'
        Write (nout, *)
        Flush (nout)
        ifail = 0
        Call nagf_file_print_matrix_real_gen('Gen', ' ', 1, n, alphar, 1, &
          'Alpha (real):', ifail)
        Call nagf_file_print_matrix_real_gen('Gen', ' ', 1, n, alphai, 1, &
          'Alpha (imag):', ifail)
        Call nagf_file_print_matrix_real_gen('Gen', ' ', 1, n, beta, 1, &
          'Beta:', ifail)
      End If
100   Continue

110   Format (1X, A, I4)
    End Program