Keyword: ヘッセンベルグ, シュール分解
概要
本サンプルは一般化上ヘッセンベルグ形の固有値の算出、または一般化シュール分解を行うFortranによるサンプルプログラムです。 本サンプルは以下に示される行列のペアに対して一般化上ヘッセンベルグ形の固有値を求めて出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン f08xef() のExampleコードです。本サンプル及びルーチンの詳細情報は f08xef のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はf08xef のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13
| このデータをダウンロード |
F08XEF Example Program Data 5 :Value of N 1.00 1.00 1.00 1.00 1.00 2.00 4.00 8.00 16.00 32.00 3.00 9.00 27.00 81.00 243.00 4.00 16.00 64.00 256.00 1024.00 5.00 25.00 125.00 625.00 3125.00 :End of matrix A 1.00 2.00 3.00 4.00 5.00 1.00 4.00 9.00 16.00 25.00 1.00 8.00 27.00 64.00 125.00 1.00 16.00 81.00 256.00 625.00 1.00 32.00 243.00 1024.00 3125.00 :End of matrix B
- １行目はタイトル行で読み飛ばされます。
- ２行目に行列A、B、Q、Zの次数(n)を指定しています。
- ３～７行目に行列Aの要素を指定しています。
- ８～１２行目に行列Bの要素を指定しています。
出力結果
（本ルーチンの詳細はf08xef のマニュアルページを参照）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
| この出力例をダウンロード |
F08XEF Example Program Results
Matrix A after balancing
1 2 3 4 5
1 1.0000 1.0000 0.1000 0.1000 0.1000
2 2.0000 4.0000 0.8000 1.6000 3.2000
3 0.3000 0.9000 0.2700 0.8100 2.4300
4 0.4000 1.6000 0.6400 2.5600 10.2400
5 0.5000 2.5000 1.2500 6.2500 31.2500
Matrix B after balancing
1 2 3 4 5
1 1.0000 2.0000 0.3000 0.4000 0.5000
2 1.0000 4.0000 0.9000 1.6000 2.5000
3 0.1000 0.8000 0.2700 0.6400 1.2500
4 0.1000 1.6000 0.8100 2.5600 6.2500
5 0.1000 3.2000 2.4300 10.2400 31.2500
Matrix A in Hessenberg form
1 2 3 4 5
1 -2.1898 -0.3181 2.0547 4.7371 -4.6249
2 -0.8395 -0.0426 1.7132 7.5194 -17.1850
3 0.0000 -0.2846 -1.0101 -7.5927 26.4499
4 0.0000 0.0000 0.0376 1.4070 -3.3643
5 0.0000 0.0000 0.0000 0.3813 -0.9937
Matrix B is triangular
1 2 3 4 5
1 -1.4248 -0.3476 2.1175 5.5813 -3.9269
2 0.0000 -0.0782 0.1189 8.0940 -15.2928
3 0.0000 0.0000 1.0021 -10.9356 26.5971
4 0.0000 0.0000 0.0000 0.5820 -0.0730
5 0.0000 0.0000 0.0000 0.0000 0.5321
Minimal required LWORK = 5
Actual value of LWORK = 30
Generalized eigenvalues
1 ( -2.437, 0.000)
2 ( 0.607, 0.795)
3 ( 0.607, -0.795)
4 ( 1.000, 0.000)
5 ( -0.410, 0.000)
- ３～８行目にバランス化された後の行列Aが出力されています。
- １１～１６行目にバランス化された後の行列Bが出力されています。
- １９～２４行目にヘッセンベルグ形の行列Aが出力されています。
- ２７～３２行目に三角行列Bが出力されています。
- ３４行目に最適なパフォーマンスを得るのに必要な最小限の作業域の配列の次元が出力されています。
- ３５行目に本関数を呼び出したプログラムの中で宣言された作業域の配列の次元の値が出力されています。
- ３８～４２行目に一般化固有値が出力されています。
ソースコード
（本ルーチンの詳細はf08xef のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
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 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
| このソースコードをダウンロード |
PROGRAM f08xefe
! F08XEF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : dgeqrf, dggbal, dgghrd, dhgeqz, dormqr, nag_wp, &
x04caf
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
INTEGER :: i, ifail, ihi, ilo, info, irows, &
jwork, lda, ldb, ldq, ldz, lwork, n
CHARACTER (1) :: compq, compz, job
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), alphai(:), alphar(:), &
b(:,:), beta(:), lscale(:), q(:,:), &
rscale(:), tau(:), work(:), z(:,:)
! .. Intrinsic Functions ..
INTRINSIC nint
! .. Executable Statements ..
WRITE (nout,*) 'F08XEF Example Program Results'
FLUSH (nout)
! Skip heading in data file
READ (nin,*)
READ (nin,*) n
ldq = 1
ldz = 1
lda = n
ldb = n
lwork = 6*n
ALLOCATE (alphai(n),alphar(n),beta(n),a(lda,n),lscale(n),q(ldq,ldq), &
rscale(n),b(ldb,n),tau(n),work(lwork),z(ldz,ldz))
! READ matrix A from data file
READ (nin,*) (a(i,1:n),i=1,n)
! READ matrix B from data file
READ (nin,*) (b(i,1:n),i=1,n)
! Balance matrix pair (A,B)
job = 'B'
! The nAG name equivalent of dggbal is f08whf
CALL dggbal(job,n,a,lda,b,ldb,ilo,ihi,lscale,rscale,work,info)
! Matrix A after balancing
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL x04caf('General',' ',n,n,a,lda,'Matrix A after balancing',ifail)
WRITE (nout,*)
FLUSH (nout)
! Matrix B after balancing
ifail = 0
CALL x04caf('General',' ',n,n,b,ldb,'Matrix B after balancing',ifail)
WRITE (nout,*)
FLUSH (nout)
! Reduce B to triangular form using QR
irows = ihi + 1 - ilo
! The nAG name equivalent of dgeqrf is f08aef
CALL dgeqrf(irows,irows,b(ilo,ilo),ldb,tau,work,lwork,info)
! Apply the orthogonal transformation to matrix A
! The nAG name equivalent of dormqr is f08agf
CALL dormqr('L','T',irows,irows,irows,b(ilo,ilo),ldb,tau,a(ilo,ilo), &
lda,work,lwork,info)
! Compute the generalized Hessenberg form of (A,B) -> (H,T)
compq = 'N'
compz = 'N'
! The nAG name equivalent of dgghrd is f08wef
CALL dgghrd(compq,compz,irows,1,irows,a(ilo,ilo),lda,b(ilo,ilo),ldb,q, &
ldq,z,ldz,info)
! Matrix A (H) in generalized Hessenberg form.
ifail = 0
CALL x04caf('General',' ',n,n,a,lda,'Matrix A in Hessenberg form', &
ifail)
WRITE (nout,*)
FLUSH (nout)
! Matrix B (T) in generalized Hessenberg form.
ifail = 0
CALL x04caf('General',' ',n,n,b,ldb,'Matrix B is triangular',ifail)
! Routine DHGEQZ
! Workspace query: jwork = -1
jwork = -1
job = 'E'
! The nAG name equivalent of dhgeqz is f08xef
CALL dhgeqz(job,compq,compz,n,ilo,ihi,a,lda,b,ldb,alphar,alphai,beta,q, &
ldq,z,ldz,work,jwork,info)
WRITE (nout,*)
WRITE (nout,99999) nint(work(1))
WRITE (nout,99998) lwork
WRITE (nout,*)
! Compute the generalized Schur form
! if the workspace lwork is adequate
IF (nint(work(1))<=lwork) THEN
! The nAG name equivalent of dhgeqz is f08xef
CALL dhgeqz(job,compq,compz,n,ilo,ihi,a,lda,b,ldb,alphar,alphai, &
beta,q,ldq,z,ldz,work,lwork,info)
! Print the generalized eigenvalues
WRITE (nout,99997)
DO i = 1, n
IF (beta(i)/=0.0E0_nag_wp) THEN
WRITE (nout,99996) i, '(', alphar(i)/beta(i), ',', &
alphai(i)/beta(i), ')'
ELSE
WRITE (nout,99994) i
END IF
END DO
ELSE
WRITE (nout,99995)
END IF
99999 FORMAT (1X,'Minimal required LWORK = ',I6)
99998 FORMAT (1X,'Actual value of LWORK = ',I6)
99997 FORMAT (1X,'Generalized eigenvalues')
99996 FORMAT (1X,I4,5X,A,F7.3,A,F7.3,A)
99995 FORMAT (1X,'Insufficient workspace allocated for call to F08XEF/DHGEQZ' &
)
99994 FORMAT (1X,I4,'Eigenvalue is infinite')
END PROGRAM f08xefe