計算ルーチン: 直交同等変換により実行列ペアの一般化実シュール分解を再構成 : (オプショナルで固有値と固有スペースの相反条件数の推定値)

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

ホーム > LAPACKサンプルプログラム目次 > 計算ルーチン > 直交同等変換により実行列ペアの一般化実シュール分解を再構成

概要

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

入力データ

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

このデータをダウンロード
DTGSEN Example Program Data
  4                    :Value of N
  4.0  1.0  1.0  2.0
  0.0  3.0  4.0  1.0
  0.0  1.0  3.0  1.0
  0.0  0.0  0.0  6.0   :End of matrix A
  2.0  1.0  1.0  3.0
  0.0  1.0  2.0  1.0
  0.0  0.0  1.0  1.0
  0.0  0.0  0.0  2.0   :End of matrix B
  1.0  0.0  0.0  0.0
  0.0  1.0  0.0  0.0
  0.0  0.0  1.0  0.0
  0.0  0.0  0.0  1.0   :End of matrix Q
  1.0  0.0  0.0  0.0
  0.0  1.0  0.0  0.0
  0.0  0.0  1.0  0.0
  0.0  0.0  0.0  1.0   :End of matrix Z
  T    F    F    T     :End of SELECT

出力結果

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

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

Warning: Floating underflow occurred
 Number of selected eigenvalues  =    2

 Selected Generalized Eigenvalues

  1 ( 2.0000E+00, 0.0000E+00)
  2 ( 3.0000E+00, 0.0000E+00)

 Norm estimate of projection onto left  eigenspace for selected cluster
   2.69E+00

 Norm estimate of projection onto right eigenspace for selected cluster
   1.50E+00

 F-norm based upper bound on Difu
   2.52E-01

 F-norm based upper bound on Difl
   2.45E-01

ソースコード

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

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


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

!     DTGSEN Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: dtgsen
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: pl, pr
      Integer :: i, ijob, info, lda, ldb, ldc, ldq, ldz, liwork, lwork, m, n
      Logical :: wantq, wantz
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
        beta(:), c(:, :), q(:, :), work(:), z(:, :)
      Real (Kind=dp) :: dif(2)
      Integer, Allocatable :: iwork(:)
      Logical, Allocatable :: select(:)
!     .. Executable Statements ..
      Write (nout, *) 'DTGSEN Example Program Results'
      Write (nout, *)
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      ldb = n
      ldc = n
      ldq = n
      ldz = n
      liwork = (n*n)/2 + 6
      lwork = n*(n+4) + 16
      Allocate (a(lda,n), alphai(n), alphar(n), b(ldb,n), beta(n), c(ldc,n), &
        q(ldq,n), work(lwork), z(ldz,n), iwork(liwork), select(n))

!     Read A, B, Q, Z and the logical array SELECT from data file

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

      Read (nin, *) select(1:n)

!     Set ijob, wantq and wantz
      ijob = 4
      wantq = .True.
      wantz = .True.

!     Reorder the Schur factors A and B and update the matrices
!     Q and Z

      Call dtgsen(ijob, wantq, wantz, select, n, a, lda, b, ldb, alphar, &
        alphai, beta, q, ldq, z, ldz, m, pl, pr, dif, work, lwork, iwork, &
        liwork, info)

      If (info>0) Then
        Write (nout, 100) info
        Write (nout, *)
        Flush (nout)
      End If

!     Print Results
      Write (nout, 130) 'Number of selected eigenvalues  = ', m
      Write (nout, *)
      Write (nout, *) 'Selected Generalized Eigenvalues'
      Write (nout, *)
      Write (nout, 120)(i, alphar(i)/beta(i), alphai(i)/beta(i), i=1, m)

      Write (nout, *)
      Write (nout, 110) 'Norm estimate of projection onto', &
        ' left  eigenspace for selected cluster', 1.0_dp/pl
      Write (nout, *)
      Write (nout, 110) 'Norm estimate of projection onto', &
        ' right eigenspace for selected cluster', 1.0_dp/pr
      Write (nout, *)
      Write (nout, 110) 'F-norm based upper bound on', ' Difu', dif(1)
      Write (nout, *)
      Write (nout, 110) 'F-norm based upper bound on', ' Difl', dif(2)

100   Format (' Reordering could not be completed. INFO = ', I3)
110   Format (1X, 2A, /, 1X, 1P, E10.2)
120   Format (1X, I2, 1X, '(', 1P, E11.4, ',', E11.4, ')')
130   Format (1X, A, I4)
    End Program


ご案内
関連情報
Privacy Policy  /  Trademarks