実一般化非対称固有値問題: 非対称行列ペア : (一般化固有値およびオプショナルで左および右固有ベクトルおよび条件数の逆数)

LAPACKサンプルソースコード : 使用ルーチン名:DGGEVX

概要

本サンプルはFortran言語によりLAPACKルーチンDGGEVXを利用するサンプルプログラムです。

行列対 $ (A,B)$の全ての固有値と右固有ベクトルを求めます。

\begin{displaymath}
A = \left(
\begin{array}{rrrr}
3.9 & 12.5 & -34.5 & -0.5 ...
...4.0 & 3.0 \\
1.0 & 3.0 & -4.0 & 4.0
\end{array} \right), \;
\end{displaymath}

条件数の推定値とそれぞれの固有値と固有ベクトルの前方誤差限界も合わせて求めます。行列対を均衡化するオプションが指定されます。

入力データ

(本ルーチンの詳細はDGGEVX のマニュアルページを参照)

このデータをダウンロード
DGGEVX 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

出力結果

(本ルーチンの詳細はDGGEVX のマニュアルページを参照)

この出力例をダウンロード
Warning: Floating underflow occurred
 DGGEVX Example Program Results


 Eigenvalue( 1) =     2.00000

 Eigenvector( 1)
     0.99606
     0.00569
     0.06261
     0.06261


 Eigenvalue( 2) = (    3.00000,   -4.00000)

 Eigenvector( 2)
 (    0.94491,   -0.00000)
 (    0.18898,    0.00000)
 (    0.11339,    0.15119)
 (    0.11339,    0.15119)


 Eigenvalue( 3) = (    3.00000,    4.00000)

 Eigenvector( 3)
 (    0.94491,    0.00000)
 (    0.18898,   -0.00000)
 (    0.11339,   -0.15119)
 (    0.11339,   -0.15119)


 Eigenvalue( 4) =     4.00000

 Eigenvector( 4)
     0.98752
     0.01097
    -0.03292
     0.15361

ソースコード

(本ルーチンの詳細はDGGEVX のマニュアルページを参照)

※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。


このソースコードをダウンロード
    Program dggevx_example

!     DGGEVX Example Program Text

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

!     .. Use Statements ..
      Use blas_interfaces, Only: dnrm2
      Use lapack_example_aux, Only: nagf_sort_realvec_rank, &
        nagf_sort_realvec_rank_rearrange
      Use lapack_interfaces, Only: dggevx
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
      Logical, Parameter :: verbose = .False.
!     .. Local Scalars ..
      Complex (Kind=dp) :: eig
      Real (Kind=dp) :: abnrm, bbnrm, eps, jswap, rcnd, scal_i, scal_r, small
      Integer :: i, ifail, ihi, ilo, info, j, k, lda, ldb, ldvr, lwork, n
      Logical :: pair
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
        beta(:), lscale(:), rconde(:), rcondv(:), rscale(:), vr(:, :), &
        vr_row(:), work(:)
      Real (Kind=dp) :: dummy(1, 1)
      Integer, Allocatable :: iwork(:)
      Logical, Allocatable :: bwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, all, cmplx, epsilon, max, maxloc, nint, sqrt, sum, &
        tiny
!     .. Executable Statements ..
      Write (nout, *) 'DGGEVX Example Program Results'
!     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), lscale(n), &
        rconde(n), rcondv(n), rscale(n), vr(ldvr,n), iwork(n+6), bwork(n))

!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call dggevx('Balance', 'No vectors (left)', 'Vectors (right)', &
        'Both reciprocal condition numbers', n, a, lda, b, ldb, alphar, &
        alphai, beta, dummy, 1, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, &
        bbnrm, rconde, rcondv, dummy, lwork, iwork, bwork, info)

!     Make sure that there is enough workspace for block size nb.
      lwork = max((nb+2*n)*n, 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 dggevx('Balance', 'No vectors (left)', 'Vectors (right)', &
        'Both reciprocal condition numbers', n, a, lda, b, ldb, alphar, &
        alphai, beta, dummy, 1, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, &
        bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)

      If (info>0) Then
        Write (nout, *)
        Write (nout, 100) 'Failure in DGGEVX. INFO =', info
      Else

!       Compute the machine precision and the safe range parameter
!       small

        eps = epsilon(1.0E0_dp)
        small = tiny(1.0E0_dp)

!       If beta(:) > eps, Order eigenvalues by ascending real parts
        If (all(abs(beta(1:n))>eps)) Then
          work(1:n) = alphar(1:n)/beta(1:n)
          ifail = 0
          Call nagf_sort_realvec_rank(work, 1, n, 'Ascending', iwork, ifail)
          Call nagf_sort_realvec_rank_rearrange(alphar, 1, n, iwork, ifail)
          Call nagf_sort_realvec_rank_rearrange(alphai, 1, n, iwork, ifail)
          Call nagf_sort_realvec_rank_rearrange(beta, 1, n, iwork, ifail)
!         Order the eigenvectors in the same way
          Allocate (vr_row(n))
          Do j = 1, n
            vr_row(1:n) = vr(j, 1:n)
            Call nagf_sort_realvec_rank_rearrange(vr_row, 1, n, iwork, ifail)
            vr(j, 1:n) = vr_row(1:n)
          End Do
          Deallocate (vr_row)
        End If

!       Print out eigenvalues and vectors and associated condition
!       number and bounds

        pair = .False.
        Do j = 1, n

!         Print out information on the j-th eigenvalue

          Write (nout, *)
          If ((abs(alphar(j))+abs(alphai(j)))*small>=abs(beta(j))) Then
            Write (nout, 110) 'Eigenvalue(', j, ')', &
              ' is numerically infinite or undetermined', 'ALPHAR(', j, &
              ') = ', alphar(j), ', ALPHAI(', j, ') = ', alphai(j), ', BETA(', &
              j, ') = ', beta(j)
          Else
            If (.Not. pair) Then
              jswap = 1.0_dp
              If (alphai(j)>0.0_dp) Then
                jswap = -jswap
              End If
            End If
            If (alphai(j)==0.0E0_dp) Then
              Write (nout, 120) 'Eigenvalue(', j, ') = ', alphar(j)/beta(j)
            Else
              eig = cmplx(alphar(j), jswap*alphai(j), kind=dp)/ &
                cmplx(beta(j), kind=dp)
              Write (nout, 130) 'Eigenvalue(', j, ') = ', eig
            End If
          End If

          If (verbose) Then
            rcnd = rconde(j)
            Write (nout, *)
            Write (nout, 140) '  Reciprocal condition number = ', rcnd
          End If

!         Print out information on the j-th eigenvector

          Write (nout, *)
          Write (nout, 150) 'Eigenvector(', j, ')'
          If (alphai(j)==0.0E0_dp) Then
!           Let largest element be positive
            work(1:n) = abs(vr(1:n,j))
            k = maxloc(work(1:n), 1)
            If (vr(k,j)<0.0_dp) Then
              vr(1:n, j) = -vr(1:n, j)/dnrm2(n, vr(1,j), 1)
            End If
            Write (nout, 160)(vr(i,j), i=1, n)
          Else
            If (pair) Then
              Write (nout, 170)(vr(i,j-1), -jswap*vr(i,j), i=1, n)
            Else
!             Let largest element be real (and positive).
              work(1:n) = vr(1:n, j)**2 + vr(1:n, j+1)**2
              k = maxloc(work(1:n), 1)
              scal_r = vr(k, j)/sqrt(work(k))/sqrt(sum(work(1:n)))
              scal_i = -vr(k, j+1)/sqrt(work(k))/sqrt(sum(work(1:n)))
              work(1:n) = vr(1:n, j)
              vr(1:n, j) = scal_r*work(1:n) - scal_i*vr(1:n, j+1)
              vr(1:n, j+1) = scal_r*vr(1:n, j+1) + scal_i*work(1:n)
              vr(k, j+1) = 0.0_dp
              Write (nout, 170)(vr(i,j), jswap*vr(i,j+1), i=1, n)
            End If
            pair = .Not. pair
          End If
          If (verbose) Then
            rcnd = rcondv(j)
            Write (nout, *)
            Write (nout, 140) '  Reciprocal condition number = ', rcnd
          End If
        End Do

      End If

100   Format (1X, A, I4)
110   Format (/, 1X, A, I2, 2A, /, 1X, 2(A,I2,A,F11.5), A, I2, A, F11.5)
120   Format (/, 1X, A, I2, A, F11.5)
130   Format (/, 1X, A, I2, A, '(', F11.5, ',', F11.5, ')')
140   Format (1X, A, 1P, E8.1)
150   Format (1X, A, I2, A)
160   Format (1X, F11.5)
170   Format (1X, '(', F11.5, ',', F11.5, ')')
    End Program


ご案内
関連情報
Privacy Policy  /  Trademarks