対称実数密行列の固有値、固有ベクトルの求解
対称実数密行列の固有値、固有ベクトルを求めるためのサンプルプログラム( tests/ABCLibTestDRSAllEigVec.f )は、以下のとおりです。
program main
include 'ABCLibDRSSED.h'
common /ABCLibpval/ nprocs, myid, ncelx, ncely, nx, ny, ierrcode
c === real*8 definition
real*8 A(0:ABCLibMAXNX-1, 0:ABCLibMAXNY-1)
real*8 eig(1:ABCLibMAXNDIVP)
real*8 X(1:ABCLibMAXN, 1:ABCLibMAXNDIVP)
c === for checking answer
real*8 dtemp
real*8 EigErr, VecErr, SysErr
c === for measuring time
real*8 t1, t2, bt, t_all
c ============================================================
c === integer definition
c === for parallel control
integer ncelx, ncely, nx, ny
integer myid, nmyidx, nmyidy
integer n, m, ms, me, nprocs
c === for error flag
integer ierrcode
c === for selecting test matrices
integer isw_mat
c ============================================================
c === Program start
c === MPI Init.
call MPI_INIT( ierr )
call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
call MPI_COMM_SIZE( MPI_COMM_WORLD, nprocs, ierr )
c =====================================================
c === Print title
if (myid .eq. 0) then
print *, ""
print *, "=== Parallel Eigensolver Check Program ==="
print *, "Parallel Eigensolver on MPI"
print *, "ABCLibDRSAllEigVec"
call ABCLibPrintVer()
print *, ""
endif
c =====================================================
c === Setting Test Condition
if (myid .eq. 0) then
write(6,*) "Problem size >"
read(5,*) n
write(6,*) "Start number of eigenvectors >"
read(5,*) ms
write(6,*) "End number of eigenvectors >"
read(5,*) me
m = me - ms + 1
c === set matrix type
isw_mat = 1
endif
c =====================================================
c == Set Parameters
call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
call MPI_BCAST(m,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
call MPI_BCAST(ms,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
call MPI_BCAST(me,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
call MPI_BCAST(isw_mat,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
c =====================================================
c === Set Parallel Control Values
if (myid .eq. 0) then
write(6,*) " ncelx="
read(5,*) ncelx
write(6,*) " ncely ="
read(5,*) ncely
endif
call MPI_BCAST(ncelx, 1、 MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
call MPI_BCAST(ncely, 1、 MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
call cid2xy2(myid, ncelx, ncely, nmyidx, nmyidy)
nx = ceiling(dfloat(n)/dfloat(ncelx))
ny = ceiling(dfloat(n)/dfloat(ncely))
c =====================================================
c === Generate Test Matrix
call MakeTestMat2(A, n, nmyidx, nmyidy, isw_mat)
c =====================================================
c === Initialize ABCLibTriRed routine
call ABCLibDRSAllEigVec_Init()
c =====================================================
c === Call Tridiagonalication Routine
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t1 = MPI_WTIME()
call ABCLibDRSAllEigVec(A, n, eig, X, ms, me)
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t2 = MPI_WTIME()
t_all = t2 - t1
call MPI_REDUCE(t_all, bt, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0,
& MPI_COMM_WORLD, ierr)
t_all = bt
c =====================================================
c === Calculate Eigenvalues & Check Its Error
if (ierrcode .eq. 0) then
if (myid .eq. 0) print *, "Normal return"
else
print *, "myid=",myid,"Abnormal return: error
code",ierrcode
endif
if (isw_mat .eq. 1) then
if (myid .eq. 0) print *、
& "Check of maximal error for eigenvalues...."
call ChkEigErr(eig, n, EigErr, ms, me)
if (myid .eq. 0) print *, "End of checking."
endif
if (myid .eq. 0) print *、
& "Check of error for eigenvectors…."
call TestVec(X, n, VecErr, m)
if (myid .eq. 0) print *, "End of checking."
if (myid .eq. 0) print *,
& "Check of error for the eigensystem…."
call TestEigSys(eig, X, n, isw_mat, SysErr, m)
if (myid .eq. 0) print *, "End of checking."
c =====================================================
c === Output Results
if (myid .eq. 0) then
print *, ""
print *, "=== Result "
print *, "================================================
"
print *, "n = ", n, " ,Number of processors
= ", nprocs
print *, "Number of eigenvalues and eigenvectors = ",
m
print *, "Start Eigenvalue No.", ms, "/ End
Eigenvalue No.",me
print *, "Test Matrix:"
if (isw_mat .eq. 1) print *, "Frank"
if (isw_mat .eq. 2) print *, "Rnd (-32768〜 32765)"
if (isw_mat .eq. 4) print *, "Glud Wilkinson W_21, d=1.0d-14"
if (isw_mat .eq. 5) print *, "Unit Matrix"
if (isw_mat .eq. 6) print *, "Wilkinson W_n+"
print *, ""
if (isw_mat .eq. 1) then
print *, "Eigenvalue Max Error = ", EigErr
endif
print *, "Eigenvector Error of ||X^T X - I|| = ",
VecErr
print *,
& "Eigensystem Error max_i||A x_i - Λ_i x_i|| = ",
SysErr
write(*, 1000) t_all
1000 format(" Total Calculation Time = ",F10.3,"
[sec]")
print *, ""
endif
c =====================================================
c === Finalize ABCLibDRSAllEig routine
call ABCLibDRSAllEigVec_Finalize()
c =====================================================
c ===== MPI finazize
call MPI_FINALIZE(ierr)
c =====================================================
c =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
stop
end
エルミート行列の全固有値、全固有ベクトルの求解
エルミート行列の全固有値、全固有ベクトルを求めるためのサンプルプログラム (tests/ABCLibTestHerAllEigVec.f)は、以下のとおりです。
program main
include 'ABCLibDRSSED.h'
common /ABCLibpval/ nprocs, myid, ncelx, ncely, nx, ny, ierrcode
c === real*8 definition
c === for tridiagonalzation
c === The AR is real part, and the AI is imaginary part of
c the Hermitian Matrix of A
real*8 AR(1:ABCLibMAXN, 1:ABCLibMAXN)
real*8 AI(1:ABCLibMAXN, 1:ABCLibMAXN)
c === The number of eigenvalues of the matrix A is 2 times
c larger than original one, because of conjugate eigenvectors.
real*8 eig(1:ABCLibMAXNDIVP*2)
real*8 X(1:ABCLibMAXN*2, 1:ABCLibMAXNDIVP*2)
c === for checking answer
real*8 dtemp
real*8 EigErr, VecErr1, VecErr2, SysErr
c === for measuring time
real*8 t1, t2, bt, t_all
c ============================================================
c === integer definition
c === for parallel control
integer ncelx, ncely, nx, ny
integer myid, nmyidx, nmyidy
integer n, nprocs
c === for error flag
integer ierrcode
c === for selecting test matrices
integer isw_mat
c ============================================================
c === Program start
c === MPI Init.
call MPI_INIT( ierr )
call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
call MPI_COMM_SIZE( MPI_COMM_WORLD, nprocs, ierr )
c =====================================================
c === Print title
if (myid .eq. 0) then
print *, ""
print *, "=== Parallel Eigensolver Check Program ==="
print *, "Parallel Eigensolver on MPI"
print *, "ABCLibHerAllEigVec"
call ABCLibPrintVer()
print *, ""
endif
c =====================================================
c === Setting Test Condition
c === set matrix type
isw_mat = 1
call MPI_BCAST(isw_mat,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
if (isw_mat .ge. 5) then
if (myid .eq. 0) then
write(6,*) "Problem size of the Hermitian
matrix >"
read(5,*) n
endif
endif
c =====================================================
c === Set Parallel Control Values
if (myid .eq. 0) then
write(6,*) " ncelx ="
read(5,*) ncelx
write(6,*) " ncely ="
read(5,*) ncely
endif
call MPI_BCAST(ncelx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
call MPI_BCAST(ncely, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
c =====================================================
c === Generate Test Matrix
call HerTestMatDim(n, isw_mat)
call HerMakeTestMat(AR, AI, n, isw_mat)
c =====================================================
c === Set Parallel Control Values
call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
nx = ceiling(dfloat(2*n)/dfloat(ncelx))
ny = ceiling(dfloat(2*n)/dfloat(ncely))
c =====================================================
c === Initialize ABCLibHerAllEig routine
call ABCLibHerAllEigVec_Init()
c =====================================================
c === Call Hermitian Eigenvalue Computation Routine
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t1 = MPI_WTIME()
call ABCLibHerAllEigVec(AR, AI, n, eig, X)
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t2 = MPI_WTIME()
t_all = t2 - t1
call MPI_REDUCE(t_all, bt, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0,
& MPI_COMM_WORLD, ierr)
t_all = bt
c =====================================================
c === Check Its Error
if (ierrcode .eq. 0) then
print *, myid, "Normal return"
else
print *, "myid=",myid,"Abnormal return: error
code",ierrcode
endif
if (isw_mat .le. 3) then
if (myid .eq. 0) then
print *,
& "Check of maximal error for eigenvalues...."
endif
call HerChkEigErr(eig, n, EigErr, isw_mat)
if (myid .eq. 0) then
if (myid .eq. 0) print *, "End of checking."
endif
endif
if (myid .eq. 0) print *,
& "Check of error for eigenvectors...."
call TestVec(X, 2*n, VecErr1, 2*n)
call TestHerVec(X, eig, n, VecErr2)
if (myid .eq. 0) print *, "End of checking."
if (myid .eq. 0) print *,
& "Check of error for the eigensystem...."
call TestHerEigSys(AR, AI, eig, X, n, isw_mat, SysErr)
if (myid .eq. 0) print *, "End of checking."
c =====================================================
c === Output Results
if (myid .eq. 0) then
print *, ""
print *, "=== Result "
print *, "================================================
"
print *, "n = ", n, " ,Number of processors
= ", nprocs
print *, "Test Matrix No:", isw_mat
print *, ""
if (isw_mat .le. 3) then
print *, "Eigenvalue Max Error = ", EigErr
endif
print *, "Eigenvector Error of Real ||X^T X - I|| = ",
VecErr1
print *, "Eigenvector Error of Her ||X^H X - I|| = ",
VecErr2
print *, "Eigensystem Error of Hermitian:"
print *, " max_i||A x_i - Λ_i x_i|| = ", SysErr
write(*, 1000) t_all
1000 format(" Total Calculation Time = ",F10.3,"
[sec]")
print *, ""
endif
c =====================================================
c === Finalize ABCLibHerAllEig routine
call ABCLibHerAllEigVec_Finalize()
c =====================================================
c ===== MPI finazize
call MPI_FINALIZE(ierr)
c =====================================================
c =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
stop
end
Gram-Schmidt法による直交化
Gram-Schmidt法による直交化をするためのサンプルプログラム (tests/ABCLibTestQRD.f) は,以下のとおりです。
program main
include 'ABCLibDRSSED.h'
common /ABCLibpval/ nprocs, myid, ncelx, ncely, nx, ny, ierrcode
c === real*8 definition
real*8 X(1:ABCLibMAXN, 1:ABCLibMAXNDIVP)
c === for checking answer
real*8 dtemp
real*8 VecErr
c === for measuring time
real*8 t1, t2, bt, t_all
c ============================================================
c === integer definition
c === for parallel control
integer ncelx, ncely, nx, ny
integer myid, nmyidx, nmyidy
integer n, nprocs
c === for error flag
integer ierrcode
c === for selecting test matrices
integer isw_mat
c ============================================================
c === Program start
c === MPI Init.
call MPI_INIT( ierr )
call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
call MPI_COMM_SIZE( MPI_COMM_WORLD, nprocs, ierr )
c =====================================================
c === Print title
if (myid .eq. 0) then
print *, ""
print *, "=== Parallel Eigensolver Check Program ==="
print *, "Parallel QR Decomposition on MPI"
print *, "ABCLibDRSAllEigVec"
call ABCLibPrintVer()
print *, ""
endif
c =====================================================
c === Setting Test Condition
if (myid .eq. 0) then
write(6,*) "problem size >"
read(5,*) n
c === set matrix type
isw_mat = 1
endif
c =====================================================
c == Set Parameters
call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
call MPI_BCAST(isw_mat,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
c =====================================================
c === Set Parallel Control Values
if (myid .eq. 0) then
write(6,*) " ncelx="
read(5,*) ncelx
write(6,*) " ncely ="
read(5,*) ncely
endif
call MPI_BCAST(ncelx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
call MPI_BCAST(ncely, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
call cid2xy2(myid, ncelx, ncely, nmyidx, nmyidy)
nx = ceiling(dfloat(n)/dfloat(ncelx))
ny = ceiling(dfloat(n)/dfloat(ncely))
c =====================================================
c === Generate Test Matrix
call MakeMatSeqCB(X, n, isw_mat)
c =====================================================
c === Initialize ABCLib_QRD routine
call ABCLibQRD_Init()
c =====================================================
c === Call Tridiagonalication Routine
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t1 = MPI_WTIME()
call ABCLib_QRD(X, n)
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t2 = MPI_WTIME()
t_all = t2 - t1
call MPI_REDUCE(t_all, bt, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0,
& MPI_COMM_WORLD, ierr)
t_all = bt
c =====================================================
c === Orthogonalized Vectors & Check Its Error
if (ierrcode .eq. 0) then
if (myid .eq. 0) print *, "Normal return"
else
print *, "myid=",myid,"Abnormal return: error
code",ierrcode
endif
if (myid .eq. 0) print *,
& "Check of error for vectors...."
call TestVec(X, n, VecErr, n)
if (myid .eq. 0) print *, "End of checking."
c =====================================================
c === Output Results
if (myid .eq. 0) then
print *, ""
print *, "=== Result "
print *, "================================================
"
print *, "n = ", n, " ,Number of processors
= ", nprocs
print *, "Test Matrix:"
if (isw_mat .eq. 1) print *, "Random"
print *, ""
print *, "Orthogonal Error of ||X^T X - I|| = ",
VecErr
write(*, 1000) t_all
1000 format(" Total Calculation Time = ",F10.3,"
[sec]")
print *, ""
endif
c =====================================================
c === Finalize ABCLib_QRD routine
call ABCLibQRD_Finalize()
c =====================================================
c ===== MPI finazize
call MPI_FINALIZE(ierr)
c =====================================================
c =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
stop
end