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

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

概要

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

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

\begin{displaymath}
A = \left(
\begin{array}{rrrr}
-3.97 - 5.04i & -4.11 + 3....
...81 - 1.59i & 3.25 + 1.33i & 1.57 - 3.44i
\end{array} \right),
\end{displaymath}

それぞれの固有値と固有ベクトルの前方誤差限界と条件数の推定値もあわせて求めます。均衡化オプションが使用されます。固有値の条件数を求めるためには左固有ベクトルを求める必要がありますがこの例で表示はされません。

入力データ

(本ルーチンの詳細はZGEEVX のマニュアルページを参照)
1
2
3
4
5
6

このデータをダウンロード
ZGEEVX Example Program Data
  4                                                            :Value of N
(-3.97, -5.04)  (-4.11,  3.70)  (-0.34,  1.01)  ( 1.29, -0.86)
( 0.34, -1.50)  ( 1.52, -0.43)  ( 1.88, -5.38)  ( 3.36,  0.65)
( 3.31, -3.85)  ( 2.50,  3.45)  ( 0.88, -1.08)  ( 0.64, -1.48)
(-1.10,  0.82)  ( 1.81, -1.59)  ( 3.25,  1.33)  ( 1.57, -3.44) :End of matrix A

出力結果

(本ルーチンの詳細はZGEEVX のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

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

 Eigenvalues

          Eigenvalue           rcond    error
  1 (-6.0004E+00,-6.9998E+00)  0.9932  3.2E-15
  2 (-5.0000E+00, 2.0060E+00)  0.9964  3.2E-15
  3 ( 7.9982E+00,-9.9637E-01)  0.9814  3.3E-15
  4 ( 3.0023E+00,-3.9998E+00)  0.9779  3.3E-15

 Eigenvectors

          Eigenvector          rcond    error

  1 ( 8.4572E-01, 0.0000E+00)  8.4011    -
    (-1.7723E-02, 3.0361E-01)
    ( 8.7521E-02, 3.1145E-01)
    (-5.6147E-02,-2.9060E-01)

  2 ( 3.8655E-01,-1.7323E-01)  8.0214    -
    ( 3.5393E-01,-4.5288E-01)
    (-6.1237E-01,-0.0000E+00)
    ( 8.5928E-02, 3.2836E-01)

  3 ( 1.7297E-01,-2.6690E-01)  5.8292    -
    (-6.9242E-01,-0.0000E+00)
    (-3.3240E-01,-4.9598E-01)
    (-2.5039E-01, 1.4655E-02)

  4 ( 3.5614E-02, 1.7822E-01)  5.8292    -
    (-1.2637E-01,-2.6663E-01)
    (-1.2933E-02, 2.9657E-01)
    (-8.8982E-01,-0.0000E+00)

 Errors below 10*machine precision are not displayed

ソースコード

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

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

!     ZGEEVX Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgeevx
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: abnrm, eps, tol
      Integer :: i, ihi, ilo, info, j, lda, ldvl, ldvr, lwork, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), vl(:, :), vr(:, :), w(:), &
        work(:)
      Complex (Kind=dp) :: dummy(1)
      Real (Kind=dp), Allocatable :: rconde(:), rcondv(:), rwork(:), scale(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: epsilon, max, nint, real
!     .. Executable Statements ..
      Write (nout, *) 'ZGEEVX Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      ldvl = n
      ldvr = n
      Allocate (a(lda,n), vl(ldvl,n), vr(ldvr,n), w(n), rconde(n), rcondv(n), &
        rwork(2*n), scale(n))

!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call zgeevx('Balance', 'Vectors (left)', 'Vectors (right)', &
        'Both reciprocal condition numbers', n, a, lda, w, vl, ldvl, vr, ldvr, &
        ilo, ihi, scale, abnrm, rconde, rcondv, dummy, lwork, rwork, info)

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

!     Read the matrix A from data file

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

!     Solve the eigenvalue problem

      Call zgeevx('Balance', 'Vectors (left)', 'Vectors (right)', &
        'Both reciprocal condition numbers', n, a, lda, w, vl, ldvl, vr, ldvr, &
        ilo, ihi, scale, abnrm, rconde, rcondv, work, lwork, rwork, info)

      If (info==0) Then

!       Compute the machine precision

        eps = epsilon(1.0E0_dp)
        tol = eps*abnrm

!       Print the eigenvalues and vectors, and associated condition
!       number and bounds

        Write (nout, *)
        Write (nout, *) 'Eigenvalues'
        Write (nout, *)
        Write (nout, *) '         Eigenvalue           rcond    error'

        Do j = 1, n

!         Print information on j-th eigenvalue

          If (rconde(j)>0.0_dp) Then
            If (tol/rconde(j)<10.0_dp*eps) Then
              Write (nout, 100) j, w(j), rconde(j), '-'
            Else
              Write (nout, 110) j, w(j), rconde(j), tol/rconde(j)
            End If
          Else
            Write (nout, 100) j, w(j), rconde(j), 'Inf'
          End If

        End Do

        Write (nout, *)
        Write (nout, *) 'Eigenvectors'
        Write (nout, *)
        Write (nout, *) '         Eigenvector          rcond    error'

        Do j = 1, n

!         Print information on j-th eigenvector

          Write (nout, *)

!         Make first real part component be positive
          If (real(vr(1,j))<0.0_dp) Then
            vr(1:n, j) = -vr(1:n, j)
          End If
          If (rcondv(j)>0.0_dp) Then
            If (tol/rcondv(j)<10.0_dp*eps) Then
              Write (nout, 100) j, vr(1, j), rcondv(j), '-'
            Else
              Write (nout, 110) j, vr(1, j), rcondv(j), tol/rcondv(j)
            End If
          Else
            Write (nout, 100) j, vr(1, j), rcondv(j), 'Inf'
          End If

          Write (nout, 120) vr(2:n, j)

        End Do
        Write (nout, *)
        Write (nout, *) 'Errors below 10*machine precision are not displayed'
      Else
        Write (nout, *)
        Write (nout, 130) 'Failure in ZGEEVX. INFO =', info
      End If

100   Format (1X, I2, 1X, '(', 1P, E11.4, ',', E11.4, ')', 1X, 0P, F7.4, 4X, &
        A)
110   Format (1X, I2, 1X, '(', 1P, E11.4, ',', E11.4, ')', 1X, 0P, F7.4, 1X, &
        1P, E8.1)
120   Format (1X, 3X, '(', 1P, E11.4, ',', E11.4, ')')
130   Format (1X, A, I4)

    End Program


ご案内
関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2025
Privacy Policy  /  Trademarks