実一般化対称固有値問題: パック形式の対称定値一般化行列 : (全ての固有値およびオプショナルで固有ベクトル)

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

ホーム > LAPACKサンプルプログラム目次 > 実一般化対称固有値問題 > パック形式の対称定値一般化行列

概要

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

一般化対称固有値問題 $ A x = \lambda B x$のすべての固有値と固有ベクトルを求めます。

\begin{displaymath}
A = \left(
\begin{array}{rrrr}
0.24 & 0.39 & 0.42 & -0.16...
...6 & 0.34 \\
-0.10 & 1.09 & 0.34 & 1.18
\end{array} \right),
\end{displaymath}

$ B$の条件数の推定値、計算された固有値と固有ベクトルの誤差限界近似値もあわせて求めます。

DSPGVDの例題プログラムは$ A B x = \lambda x$の形式の一般化対称固有値問題をとく方法を示します。

入力データ

(本ルーチンの詳細はDSPGV のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13

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

  4                         :Value of N

  0.24   0.39   0.42  -0.16
        -0.11   0.79   0.63
               -0.25   0.48
                      -0.03 :End of matrix A

  4.16  -3.12   0.56  -0.10
         5.03  -0.83   1.09
                0.76   0.34
                       1.18 :End of matrix B

出力結果

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

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

 Eigenvalues
       -2.2254    -0.4548     0.1001     1.1270

 Estimate of reciprocal condition number for B
        5.8E-03

 Error estimates for the eigenvalues
        9.3E-14    2.5E-14    1.1E-14    5.1E-14

ソースコード

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

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

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

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

!     DSPGV Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: dlansp, dspgv, dtpcon
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
      Character (1), Parameter :: uplo = 'U'
!     .. Local Scalars ..
      Real (Kind=dp) :: anorm, bnorm, eps, rcond, rcondb, t1, t2
      Integer :: i, info, j, n
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: ap(:), bp(:), eerbnd(:), w(:), work(:)
      Real (Kind=dp) :: dummy(1, 1)
      Integer, Allocatable :: iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, epsilon
!     .. Executable Statements ..
      Write (nout, *) 'DSPGV Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n

      Allocate (ap((n*(n+1))/2), bp((n*(n+1))/2), eerbnd(n), w(n), work(3*n), &
        iwork(n))

!     Read the upper or lower triangular parts of the matrices A and
!     B from data file

      If (uplo=='U') Then
        Read (nin, *)((ap(i+(j*(j-1))/2),j=i,n), i=1, n)
        Read (nin, *)((bp(i+(j*(j-1))/2),j=i,n), i=1, n)
      Else If (uplo=='L') Then
        Read (nin, *)((ap(i+((2*n-j)*(j-1))/2),j=1,i), i=1, n)
        Read (nin, *)((bp(i+((2*n-j)*(j-1))/2),j=1,i), i=1, n)
      End If

!     Compute the one-norms of the symmetric matrices A and B

      anorm = dlansp('One norm', uplo, n, ap, work)
      bnorm = dlansp('One norm', uplo, n, bp, work)

!     Solve the generalized symmetric eigenvalue problem
!     A*x = lambda*B*x (itype = 1)

      Call dspgv(1, 'No vectors', uplo, n, ap, bp, w, dummy, 1, work, info)

      If (info==0) Then

!       Print solution

        Write (nout, *) 'Eigenvalues'
        Write (nout, 100) w(1:n)

!       Call DTPCON to estimate the reciprocal condition
!       number of the Cholesky factor of B.  Note that:
!       cond(B) = 1/rcond**2

        Call dtpcon('One norm', uplo, 'Non-unit', n, bp, rcond, work, iwork, &
          info)

!       Print the reciprocal condition number of B

        rcondb = rcond**2
        Write (nout, *)
        Write (nout, *) 'Estimate of reciprocal condition number for B'
        Write (nout, 110) rcondb

!       Get the machine precision, eps, and if rcondb is not less
!       than eps**2, compute error estimates for the eigenvalues

        eps = epsilon(1.0E0_dp)
        If (rcond>=eps) Then
          t1 = eps/rcondb
          t2 = anorm/bnorm
          Do i = 1, n
            eerbnd(i) = t1*(t2+abs(w(i)))
          End Do

!         Print the approximate error bounds for the eigenvalues

          Write (nout, *)
          Write (nout, *) 'Error estimates for the eigenvalues'
          Write (nout, 110) eerbnd(1:n)
        Else
          Write (nout, *)
          Write (nout, *) 'B is very ill-conditioned, error ', &
            'estimates have not been computed'
        End If
      Else If (info>n .And. info<=2*n) Then
        i = info - n
        Write (nout, 120) 'The leading minor of order ', i, &
          ' of B is not positive definite'
      Else
        Write (nout, 130) 'Failure in DSPGV. INFO =', info
      End If

100   Format (3X, (6F11.4))
110   Format (4X, 1P, 6E11.1)
120   Format (1X, A, I4, A)
130   Format (1X, A, I4)
    End Program


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