対称実数密行列の固有値、固有ベクトルの求解


対称実数密行列の固有値、固有ベクトルを求めるためのサンプルプログラム( 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