計算機ソフトウェアを利用した固有値と固有ベクトルの計算の練習
[1] GNU Octave (https://www.gnu.org/software/octave/)
[2] R (https://www.r-project.org/)
[3] Fortran (http://www.nag-j.co.jp/fortran/)とLAPACK (http://www.netlib.org/lapack/)
どちらかというと[1], [2]は練習やすぐ結果を見たいときに使うと良い。
[3]は大規模な行列や似たような処理を多数回行うときに使うと良い。
[1] GNU Octave
$ octave
GNU Octave, version 3.0.5
octave:1> A=[1.0, -0.3; -0.7, 0.6]
A =
1.00000 -0.30000
-0.70000 0.60000
octave:2> eig(A)
ans =
1.30000
0.30000
固有ベクトルと固有値を求める
octave:3> [V, LAMBDA] = eig (A)
V =
0.70711 0.39392
-0.70711 0.91915
１列目が１個目の固有ベクトル, 2列目が2個目の固有ベクトル
LAMBDA =
1.30000 0.00000
0.00000 0.30000
(1,1)成分は1個目の固有ベクトルに対応する固有値
(2,2)成分は2個目の固有ベクトルに対応する固有値
固有ベクトルは方向を変えないことの確認
octave:4> A*V
ans =
0.91924 0.11818
-0.91924 0.27574
octave:5> 0.70711*1.3
ans = 0.91924
octave:6> 0.39392*0.3
ans = 0.11818
octave:7> 0.91915*0.3
ans = 0.27575
[2] R
$ R
R version 2.10.0 (2009-10-26)
> A<-rbind(c(1,-0.3),c(-0.7,0.6))
> A
[,1] [,2]
[1,] 1.0 -0.3
[2,] -0.7 0.6
> z<-eigen(A)
> z$values
[1] 1.3 0.3
> z$vectors[,1]
[1] 0.7071068 -0.7071068
> z$vectors[,2]
[1] 0.3939193 0.9191450
固有ベクトルは方向を変えないことの確認
> A%*%z$vectors
[,1] [,2]
[1,] 0.9192388 0.1181758
[2,] -0.9192388 0.2757435
> z$values[1]*z$vectors[,1]
[1] 0.9192388 -0.9192388
> z$values[2]*z$vectors[,2]
[1] 0.1181758 0.2757435
CTL+D
Save workspace image? [y/n/c]: n
[3] Fortran with Intel Math Kernel
[2016年 12月 16日 金曜日 10:38:40 JST]
[~/TEACHING/TOOLS/Linear.Algebra/Eigin]
[am@aofd165]
$ ifort -V
インテル(R) インテル(R) 64 Fortran コンパイラー XE (インテル(R) 64 対応アプリケーション用)、バージョン 13.1.1.163 ビルド 20130313
[2016年 12月 17日 土曜日 00:44:04 JST]
[~/TEACHING/TOOLS/Linear.Algebra/Eigin]
[am@aofd165]
$ eigin.sh
Compiling eigin.f90 ...
Done Compile.
eigin.exe is running ...
DGEEV Example Program Results
Eigenvalue( 1) = 1.3000E+00
Eigenvector( 1)
0.7071E+00
-0.7071E+00
Eigenvalue( 2) = 3.0000E-01
Eigenvector( 2)
0.3939E+00
0.9191E+00
Done eigin.exe
ソースコード
$ srcdump.sh eigin.sh
------------------------------
List of the following files:
------------------------------
eigin.sh
------------------------------
Machine info
------------------------------
aofd165.bio.mie-u.ac.jp
/work1/am/TEACHING/TOOLS/Linear.Algebra/Eigin
Sat Dec 17 00:45:01 JST 2016
======================
eigin.sh
======================
#!/bin/bash
# Description:
src=$(basename $0 .sh).f90
exe=$(basename $0 .sh).exe
data_file=$(basename $0 .sh)_test.txt
f90=ifort
opt="-CB -traceback -fpe0 -mkl"
temp=$(basename $src .f90)
prog_name=$(echo $temp| sed -e 's/\./_/g')
cat <<EOF > ${data_file}
DGEEV Example Program Data
2 :Value of N
1.00 -0.30
-0.70 0.60 :End of matrix A
EOF
cat <<EOF>$src
program ${prog_name}
integer,parameter::nin=10, nout=6
integer,parameter::nb=64, nmax=10
integer,parameter::lda=nmax, ldvr=nmax, lwork=(2+nb)*nmax
complex (kind(0D0)) eig
integer i, info, j, lwkopt, n
double Precision::a(lda, nmax), dummy(1, 1), vr(ldvr, nmax), &
wi(nmax), work(lwork), wr(nmax)
external dgeev
intrinsic cmplx
write (nout, *) 'DGEEV Example Program Results'
! Skip heading in data file
open(10,file="${data_file}",action="read")
read (nin, *)
read (nin, *) n
if (n<=nmax) Then
!
! Read the matrix A from data file
!
read (nin, *)((a(i,j),j=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)
lwkopt = work(1)
!
if (info==0) Then
do j = 1, n
write (nout, *)
if (wi(j)==0.0D0) Then
write (nout,'(1X, A, I2, A, 1P, E11.4)') &
'Eigenvalue(', j, ') = ', wr(j)
else
eig = cmplx(wr(j), wi(j), kind(0D0))
write (nout, '(A,i5,A,E13.5)') 'Eigenvalue(', j, ')', eig
end If
write (nout, *)
write (nout, '(1X, A, I2, A)') 'Eigenvector(', j, ')'
if (wi(j)==0.0D0) Then
write (nout, '(1X, E11.4)') (vr(i,j), i=1, n)
else If (wi(j)>0.0D0) then
write (nout, '(E11.4, E11.4)') (vr(i,j), vr(i,j+1), i=1, n)
else
write (nout, '(E11.4, E11.4)' ) &
(vr(i,j-1), -vr(i,j), i=1, n)
end If
end Do
else
write (nout, *)
write (nout, '(1X, A, I4)') 'Failure in DGEEV. INFO = ', info
end If
!
! Print workspace information
!
if (lwork<lwkopt) Then
Write (nout, *)
Write (nout, '(1X, A, I5, /1X, A, I5)') &
'Optimum workspace required = ', lwkopt, &
'Workspace provided = ', lwork
End If
Else
Write (nout, *)
Write (nout, *) 'NMAX too small'
End If
Stop
!
stop
end program ${prog_name}
EOF
echo
echo Compiling ${src} ...
${f90} ${opt} ${src} -o ${exe}
if [ $? -ne 0 ]; then
echo
echo Compile error!
echo
echo Terminated.
echo
exit 1
fi
echo "Done Compile."
echo
echo ${exe} is running ...
echo
${exe}
if [ $? -ne 0 ]; then
echo
echo ERROR in $exe: Runtime error!
echo
echo Terminated.
echo
exit 1
fi
echo
echo "Done ${exe}"
echo
----------------------
End of eigin.sh
----------------------
EOF解析の場合、対称行列の固有値・固有ベクトルを求めることになる。その場合、下記の例のように、DSYEVを使ったほうが良い。
$ eigin.sym.sh
Compiling eigin.sym.f90 ...
Done Compile.
eigin.sym.exe is running ...
DSYEV Example Program Results
Eigenvalues
-2.0531 -0.5146 -0.2943 12.8621
Eigenvectors
0.7003 -0.5144 0.2767 0.4103
0.3592 0.4851 -0.6634 0.4422
-0.1569 0.5420 0.6504 0.5085
-0.5965 -0.4543 -0.2457 0.6144
Error estimate for the eigenvalues
1.4E-15
Error estimates for the eigenvectors
9.3E-16 6.5E-15 6.5E-15 1.1E-16
Check results:
Input matrix, A:
1.0000 2.0000 3.0000 4.0000
2.0000 2.0000 3.0000 4.0000
3.0000 3.0000 3.0000 4.0000
4.0000 4.0000 4.0000 4.0000
A*x lambda*x
-1.4379 -1.4379
-0.7375 -0.7375
0.3220 0.3220
1.2248 1.2248
0.2647 0.2647
-0.2497 -0.2497
-0.2789 -0.2789
0.2338 0.2338
-0.0814 -0.0814
0.1952 0.1952
-0.1914 -0.1914
0.0723 0.0723
5.2779 5.2779
5.6882 5.6882
6.5408 6.5408
7.9019 7.9019
Done eigin.sym.exe
$ srcdump.sh eigin.sym.sh
------------------------------
List of the following files:
------------------------------
eigin.sym.sh
------------------------------
Machine info
------------------------------
aofd165.bio.mie-u.ac.jp
/work1/am/TEACHING/TOOLS/Linear.Algebra/Eigin.Symmetric
Mon Dec 19 14:54:33 JST 2016
======================
eigin.sym.sh
======================
#!/bin/bash
# Description:
src=$(basename $0 .sh).f90
exe=$(basename $0 .sh).exe
data_file=$(basename $0 .sh)_test.txt
f90=ifort
opt="-CB -traceback -fpe0 -mkl"
temp=$(basename $src .f90)
prog_name=$(echo $temp| sed -e 's/\./_/g')
cat <<EOF > ${data_file}
DSYEV Example Program Data
4 !Value of N
1.0 2.0 3.0 4.0
2.0 3.0 4.0
3.0 4.0
4.0 !End of matrix A
EOF
cat <<EOF>$src
program ${prog_name}
integer,parameter::nin=5, nout=6
integer,parameter::nb=64
double precision::eerrbd, eps
integer::i, ifail, info, j, lwkopt, n
double precision,allocatable::a(:,:), a_in(:,:), rcondz(:), w(:), &
work(:), zerrbd(:)
double precision:: dlamch
external::ddisna, dsyev
! For checking
double precision,parameter::al=1.0d0, beta=0.0d0 ! For dgemv & dgemm
double precision,allocatable:: c(:,:),x(:),y(:)
write (nout, *) 'DSYEV Example Program Results'
write (nout, *)
!
! Read in input data
!
read (nin, *)
read (nin, *) n
lda=n
lwork=(nb+2)*n
allocate(a(lda,n), a_in(lda,n), rcondz(n), w(n), work(lwork), zerrbd(n))
allocate(c(lda,n),x(lda),y(lda))
read (nin, *)((a(i,j),j=i,n), i=1, n)
do i=1,n
do j=1,n
if(j >= i)then
a_in(i,j)=a(i,j)
else
a_in(i,j)=a(j,i)
end if
end do
end do
!
! Compute eigenvalues and eigenvectors
!
call dsyev('Vectors', 'Upper', n, a, lda, w, work, lwork, info)
lwkopt = work(1)
if (info/=0)then
print *,'Error in DSYEV. INFO =', info
stop
endif
write(nout, *) 'Eigenvalues'
write(nout, '(3X, (8F8.4))')(w(j), j=1, n)
ifail = 0
write(nout,*)
write(nout,*) 'Eigenvectors'
do i=1,n
write(nout,'(100f8.4)')(a(i,j),j=1,n)
end do !i
!
! Error estimate
!
eps = dlamch('Eps')
eerrbd = eps*max(abs(w(1)), abs(w(n)))
call ddisna('Eigenvectors', n, n, w, rcondz, info)
do i = 1, n
zerrbd(i) = eerrbd/rcondz(i)
end do
write (nout, *)
write (nout, *) 'Error estimate for the eigenvalues'
write (nout, '(4X, 1P, 6E11.1)') eerrbd
write (nout, *)
write (nout, *) 'Error estimates for the eigenvectors'
write (nout, '(4X, 1P, 6E11.1)')(zerrbd(i), i=1, n)
if (lwork<lwkopt) Then
write (nout, *)
write (nout, '(1X, A, I5, /1X, A, I5)') &
'Optimum workspace required = ', lwkopt, &
'Workspace provided = ', lwork
end if
!
! Check results
!
print *
print *,'Check results:'
c(:,:)=0.0d0
call dgemm('n','n', n, n, n, al, a_in, lda, a, lda, beta, c, lda)
print *,'Input matrix, A:'
do i=1,n
write(nout,'(100f8.4)')(a_in(i,j),j=1,n)
end do !i
print *
print '(A)',' A*x lambda*x'
do j=1,n
y(:)=w(j)*a(:,j)
do i=1,n
print '(f8.4,3x,f8.4)', c(i,j), y(i)
end do !i
print *
end do !j
stop
end program ${prog_name}
EOF
echo
echo Compiling ${src} ...
${f90} ${opt} ${src} -o ${exe}
if [ $? -ne 0 ]; then
echo
echo Compile error!
echo
echo Terminated.
echo
exit 1
fi
echo "Done Compile."
echo
echo ${exe} is running ...
echo
${exe} < ${data_file}
if [ $? -ne 0 ]; then
echo
echo ERROR in $exe: Runtime error!
echo
echo Terminated.
echo
exit 1
fi
echo
echo "Done ${exe}"
echo
----------------------
End of eigin.sym.sh
----------------------