実非対称固有値問題: 非対称行列 : (全ての固有値およびオプショナルで左および右固有ベクトル)

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

概要

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

行列のすべての固有値と右固有ベクトルを求めます。

\begin{displaymath}
A = \left(
\begin{array}{rrrr}
0.35 & 0.45 & -0.14 & -0.1...
... & 0.17 \\
0.25 & -0.32 & -0.13 & 0.11
\end{array} \right).
\end{displaymath}

入力データ

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

このデータをダウンロード
DGEEV Example Program Data

  4                         :Value of N

  0.35   0.45  -0.14  -0.17
  0.09   0.07  -0.54   0.35
 -0.44  -0.33  -0.03   0.17
  0.25  -0.32  -0.13   0.11 :End of matrix A

出力結果

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

この出力例をダウンロード
 DGEEV Example Program Results

 Eigenvalue( 1) =  7.9948E-01

 Eigenvector( 1)
 -6.5509E-01
 -5.2363E-01
  5.3622E-01
 -9.5607E-02

 Eigenvalue( 2) = (-9.9412E-02, 4.0079E-01)

 Eigenvector( 2)
 (-1.9330E-01, 2.5463E-01)
 ( 2.5186E-01,-5.2240E-01)
 ( 9.7182E-02,-3.0838E-01)
 ( 6.7595E-01, 0.0000E+00)

 Eigenvalue( 3) = (-9.9412E-02,-4.0079E-01)

 Eigenvector( 3)
 (-1.9330E-01,-2.5463E-01)
 ( 2.5186E-01, 5.2240E-01)
 ( 9.7182E-02, 3.0838E-01)
 ( 6.7595E-01,-0.0000E+00)

 Eigenvalue( 4) = -1.0066E-01

 Eigenvector( 4)
  1.2533E-01
  3.3202E-01
  5.9384E-01
  7.2209E-01

ソースコード

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

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


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

!     DGEEV Example Program Text
!     Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com

!     .. Use Statements ..
      Use lapack_interfaces, Only: dgeev
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=dp), Parameter :: zero = 0.0_dp
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Complex (Kind=dp) :: eig
      Real (Kind=dp) :: alpha, beta, scale
      Integer :: i, info, j, k, lda, ldvr, lwork, n
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), vr(:, :), wi(:), work(:), wr(:)
      Real (Kind=dp) :: dummy(1, 1)
!     .. Intrinsic Procedures ..
      Intrinsic :: cmplx, max, maxloc, nint, sqrt
!     .. Executable Statements ..
      Write (nout, *) 'DGEEV Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      ldvr = n
      Allocate (a(lda,n), vr(ldvr,n), wi(n), wr(n))

!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call dgeev('No left vectors', 'Vectors (right)', n, a, lda, wr, wi, &
        dummy, 1, vr, ldvr, dummy, lwork, info)

!     Make sure that there is enough workspace for block size nb.
      lwork = max((nb+2)*n, nint(dummy(1,1)))
      Allocate (work(lwork))

!     Read the matrix A from data file

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

!     Compute the eigenvalues and right eigenvectors of A

      Call dgeev('No left vectors', 'Vectors (right)', n, a, lda, wr, wi, &
        dummy, 1, vr, ldvr, work, lwork, info)

      If (info==0) Then

!       Print solution

        Do j = 1, n
          Write (nout, *)
          If (wi(j)==zero) Then
            Write (nout, 100) 'Eigenvalue(', j, ') = ', wr(j)
          Else
            eig = cmplx(wr(j), wi(j), kind=dp)
            Write (nout, 110) 'Eigenvalue(', j, ') = ', eig
          End If
          Write (nout, *)
          Write (nout, 120) 'Eigenvector(', j, ')'
          If (wi(j)==zero) Then
!           Scale by making largest element positive
            k = maxloc(vr(1:n,j), 1)
            If (vr(k,j)<zero) Then
              vr(1:n, j) = -vr(1:n, j)
            End If
            Write (nout, 130) vr(1:n, j)
          Else If (wi(j)>0.0E0_dp) Then
!           Scale by making largest element real and positive
            work(1:n) = vr(1:n, j)**2 + vr(1:n, j+1)**2
            k = maxloc(work(1:n), 1)
            scale = sqrt(work(k))
            work(1:n) = vr(1:n, j)
            alpha = vr(k, j)/scale
            beta = vr(k, j+1)/scale
            vr(1:n, j) = alpha*work(1:n) + beta*vr(1:n, j+1)
            vr(1:n, j+1) = alpha*vr(1:n, j+1) - beta*work(1:n)
            Write (nout, 140)(vr(i,j), vr(i,j+1), i=1, n)
          Else
            Write (nout, 140)(vr(i,j-1), -vr(i,j), i=1, n)
          End If
        End Do
      Else
        Write (nout, *)
        Write (nout, 150) 'Failure in DGEEV.  INFO = ', info
      End If

100   Format (1X, A, I2, A, 1P, E11.4)
110   Format (1X, A, I2, A, '(', 1P, E11.4, ',', 1P, E11.4, ')')
120   Format (1X, A, I2, A)
130   Format (1X, 1P, E11.4)
140   Format (1X, '(', 1P, E11.4, ',', 1P, E11.4, ')')
150   Format (1X, A, I4)
    End Program


ご案内
関連情報
Privacy Policy  /  Trademarks