Commit 8a9c9df1 authored by Pavel Kus's avatar Pavel Kus

Introducing analytical test

Introducing new test in which matrix and its eigendecomposition is
known and thus can be easily created and checked directly, without the
need to use scalapack or any other communication (apart from reducing
error).

The test is based on the fact, that if L_A and S_A are eigenvalues and
eigenvectors of matrix A, respectively, and L_B and S_B eigenvalues and
eigenvectors of B, then kron(L_A, L_B) and kron (S_A, S_B) are
eigenvalues and eigenvectors of kron(A, B).
Since it is easy to know exact eigendecomposition of a small matrix (e.g.
2x2), and kron operator has very simple structure, we can construct
arbitrarily large matrix and its eigendecomposition. We only have to
select small matrices such that the resulting matrix has unique and
ordered eigenvalues, so that the checking of the result is than easy.
Each element of matrix, eigenvector matrix and eigenvalue vector can
be quickly computed independently, just using its global coordinates.

The test is currently limited to matrices of size 2^n, but by
storing eigendecompositions of more small matrices (e.g. 3x3 and 5x5) we
could construct any matrix of size 2^n*3^m*5^o, which would probably be
sufficient, since most often used sizes (150, 1000, 5000, 2000, 60000)
are of this form.
parent dde98bdb
......@@ -459,6 +459,7 @@ libelpatest@SUFFIX@_la_SOURCES = \
test/shared/test_blacs_infrastructure.F90 \
test/shared/test_prepare_matrix.F90 \
test/shared/test_prepare_matrix_template.X90 \
test/shared/test_analytic.F90 \
test/shared/test_output_type.F90
if HAVE_REDIRECT
......
......@@ -18,6 +18,10 @@ gpu_flag = {
0 : "-DTEST_GPU=0",
1 : "-DTEST_GPU=1",
}
matrix_flag = {
"random" : "-DTEST_MATRIX_RANDOM",
"analytic" : "-DTEST_MATRIX_ANALYTIC",
}
test_type_flag = {
"eigenvectors" : "-D__EIGENVECTORS",
......@@ -25,12 +29,17 @@ test_type_flag = {
"solve_tridiagonal" : "-D__SOLVE_TRIDIAGONAL",
}
for g, t, p, d, s in product(sorted(gpu_flag.keys()),
for m, g, t, p, d, s in product(sorted(matrix_flag.keys()),
sorted(gpu_flag.keys()),
sorted(test_type_flag.keys()),
sorted(prec_flag.keys()),
sorted(domain_flag.keys()),
sorted(solver_flag.keys())):
#todo: decide what tests we actually want
if(m == "analytic" and (g == 1 or t != "eigenvectors" or p == "single" or d == "complex")):
continue
for kernel in ["all_kernels", "default_kernel"] if s == "2stage" else ["nokernel"]:
endifs = 0
extra_flags = []
......@@ -52,7 +61,7 @@ for g, t, p, d, s in product(sorted(gpu_flag.keys()),
raise Exception("Oh no!")
endifs += 1
name = "test_{0}_{1}_{2}_{3}{4}{5}".format(d, p, t, s, "" if kernel == "nokernel" else "_" + kernel, "_gpu" if g else "")
name = "test_{0}_{1}_{2}_{3}{4}{5}{6}".format(d, p, t, s, "" if kernel == "nokernel" else "_" + kernel, "_gpu" if g else "", "_analytic" if m == "analytic" else "")
print("noinst_PROGRAMS += " + name)
print("check_SCRIPTS += " + name + ".sh")
print(name + "_SOURCES = test/Fortran/test.F90")
......@@ -63,6 +72,7 @@ for g, t, p, d, s in product(sorted(gpu_flag.keys()),
prec_flag[p],
test_type_flag[t],
solver_flag[s],
gpu_flag[g]] + extra_flags))
gpu_flag[g],
matrix_flag[m]] + extra_flags))
print("endif\n" * endifs)
......@@ -110,6 +110,7 @@ program test
use test_read_input_parameters
use test_blacs_infrastructure
use test_check_correctness
use test_analytic
implicit none
......@@ -124,7 +125,7 @@ program test
integer :: mpierr
! blacs
integer :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol
integer :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol, adjusted_na
! The Matrix
MATRIX_TYPE, allocatable :: a(:,:), as(:,:)
......@@ -157,6 +158,7 @@ program test
stop 77
#endif
call read_input_parameters_traditional(na, nev, nblk, write_to_file)
call setup_mpi(myid, nprocs)
do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
......@@ -165,6 +167,19 @@ program test
np_rows = nprocs/np_cols
#ifdef TEST_MATRIX_ANALYTIC
adjusted_na = 1
do while (adjusted_na < na)
adjusted_na = adjusted_na * 2
end do
if (adjusted_na > na) then
na = adjusted_na
if(myid == 0) then
print *, 'At the moment, analytic test works for sizes of matrix of powers of two only. na changed to ', na
end if
end if
#endif
if (myid == 0) then
print '((a,i0))', 'Matrix size: ', na
print '((a,i0))', 'Num eigenvectors: ', nev
......@@ -180,7 +195,10 @@ program test
call set_up_blacs_descriptor(na, nblk, my_prow, my_pcol, np_rows, np_cols, &
na_rows, na_cols, sc_desc, my_blacs_ctxt, info)
allocate(a (na_rows,na_cols), as(na_rows,na_cols))
allocate(a (na_rows,na_cols))
#ifdef TEST_MATRIX_RANDOM
allocate(as(na_rows,na_cols))
#endif
allocate(z (na_rows,na_cols))
allocate(ev(na))
......@@ -195,8 +213,12 @@ program test
ev(:) = 0.0
#ifdef __EIGENVECTORS
#ifdef TEST_MATRIX_ANALYTIC
call prepare_matrix_analytic(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
#else
call prepare_matrix(na, myid, sc_desc, a, z, as)
#endif
#endif
#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
! set up simple toeplitz matrix
......@@ -352,7 +374,11 @@ program test
endif
#ifdef __EIGENVECTORS
#ifdef TEST_MATRIX_ANALYTIC
status = check_correctness_analytic(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
#else
status = check_correctness(na, nev, as, z, ev, sc_desc, myid)
#endif
if (status /= 0) then
if (myid == 0) print *, "Result incorrect!"
call exit(status)
......@@ -424,7 +450,11 @@ program test
#endif
#ifdef TEST_ALL_KERNELS
#ifdef TEST_MATRIX_ANALYTIC
call prepare_matrix_analytic(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
#else
a(:,:) = as(:,:)
#endif
#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
d = ds
sd = sds
......@@ -436,7 +466,9 @@ program test
call elpa_uninit()
deallocate(a)
#ifdef TEST_MATRIX_RANDOM
deallocate(as)
#endif
deallocate(z)
deallocate(ev)
......
!#include "assert.h" ! why complains?
module test_analytic
use test_util
interface prepare_matrix_analytic
module procedure prepare_matrix_analytic_real_double
end interface
interface check_correctness_analytic
module procedure check_correctness_analytic_real_double
end interface
contains
subroutine prepare_matrix_analytic_real_double (na, a, nblk, myid, np_rows, &
np_cols, my_prow, my_pcol)
use elpa_utilities
implicit none
integer(kind=ik), intent(in) :: na, nblk, myid, np_rows, np_cols, my_prow, my_pcol
real(kind=rk8), intent(inout) :: a(:,:)
integer(kind=ik) :: globI, globJ, locI, locJ, levels
if(.not. decompose(na, levels)) then
print *, "can not decomopse matrix size"
stop 1
end if
do globI = 1, na
do globJ = 1, na
if(map_global_array_index_to_local_index(globI, globJ, locI, locJ, &
nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(locI, locJ) = analytic_matrix(levels, globI, globJ)
end if
end do
end do
end subroutine
function check_correctness_analytic_real_double (na, nev, ev, z, nblk, myid, np_rows, &
np_cols, my_prow, my_pcol) result(status)
use elpa_utilities
implicit none
integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol
integer(kind=ik) :: status
real(kind=rk8), intent(inout) :: z(:,:)
real(kind=rk8), intent(inout) :: ev(:)
integer(kind=ik) :: globI, globJ, locI, locJ, levels
real(kind=rk8) :: diff, max_z_diff, max_ev_diff, glob_max_z_diff, max_curr_z_diff_minus, max_curr_z_diff_plus
real(kind=rk8) :: computed, expected
if(.not. decompose(na, levels)) then
print *, "can not decomopse matrix size"
stop 1
end if
max_z_diff = 0.0_rk8
max_ev_diff = 0.0_rk8
do globJ = 1, nev
diff = abs(ev(globJ) - analytic_eigenvalues(levels, globJ))
max_ev_diff = max(diff, max_ev_diff)
! calculated eigenvector can be in opposite direction
max_curr_z_diff_minus = 0.0_rk8
max_curr_z_diff_plus = 0.0_rk8
do globI = 1, na
if(map_global_array_index_to_local_index(globI, globJ, locI, locJ, &
nblk, np_rows, np_cols, my_prow, my_pcol)) then
computed = z(locI, locJ)
expected = analytic_eigenvectors(levels, globI, globJ)
max_curr_z_diff_minus = max(abs(computed - expected), max_curr_z_diff_minus)
max_curr_z_diff_plus = max(abs(computed + expected), max_curr_z_diff_plus)
end if
end do
! we have max difference of one of the eigenvectors, update global
max_z_diff = max(max_z_diff, min(max_curr_z_diff_minus, max_curr_z_diff_plus))
end do
#ifdef WITH_MPI
call mpi_allreduce(max_z_diff, glob_max_z_diff, 1, MPI_REAL8, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else
glob_max_z_diff = max_z_diff
#endif
if(myid == 0) print *, 'Maximal error in eigenvalues :', max_ev_diff
if(myid == 0) print *, 'Maximal error in eigenvectors :', glob_max_z_diff
status = 0
if (max_ev_diff .gt. 5e-14_rk8 .or. max_ev_diff .eq. 0.0_rk8) status = 1
if (glob_max_z_diff .gt. 6e-13_rk8 .or. glob_max_z_diff .eq. 0.0_rk8) status = 1
end function
function decompose(num, decomposition) result(possible)
implicit none
integer(kind=ik), intent(in) :: num
integer(kind=ik), intent(out) :: decomposition
logical :: possible
integer(kind=ik) :: reminder
decomposition = 0
possible = .true.
reminder = num
do while (reminder > 1)
if (MOD(reminder, 2) == 0) then
decomposition = decomposition + 1
reminder = reminder / 2
else
possible = .false.
end if
end do
end function
#define ANALYTIC_MATRIX 0
#define ANALYTIC_EIGENVECTORS 1
#define ANALYTIC_EIGENVALUES 2
function analytic_matrix(levels, i, j) result(element)
implicit none
integer(kind=ik), intent(in) :: levels, i, j
real(kind=rk8) :: element
element = analytic(levels, i, j, ANALYTIC_MATRIX)
end function
function analytic_eigenvectors(levels, i, j) result(element)
implicit none
integer(kind=ik), intent(in) :: levels, i, j
real(kind=rk8) :: element
element = analytic(levels, i, j, ANALYTIC_EIGENVECTORS)
end function
function analytic_eigenvalues(levels, i) result(element)
implicit none
integer(kind=ik), intent(in) :: levels, i
real(kind=rk8) :: element
element = analytic(levels, i, i, ANALYTIC_EIGENVALUES)
end function
function analytic(n, i, j, what) result(element)
implicit none
integer(kind=ik), intent(in) :: n, i, j, what
real(kind=rk8) :: element, am
real(kind=rk8) :: a, s, c, mat(2,2)
integer(kind=ik) :: ii, jj, m
! assert(i < 2**n)
! assert(j < 2**n)
! assert(i >= 0)
! assert(j >= 0)
! go to zero-based indexing
ii = i - 1
jj = j - 1
a = get_a(n)
s = 0.5_rk8
c = sqrt(3.0_rk8)/2.0_rk8
element = 1.0_rk8
do m = 1, n
am = a**(2**(m-1))
if(what == ANALYTIC_MATRIX) then
mat = reshape((/ c*c + am * s*s, (1.0_rk8-am) * s*c, &
(1.0_rk8-am) * s*c, s*s + am * c*c /), &
(/2, 2/))
else if(what == ANALYTIC_EIGENVECTORS) then
mat = reshape((/ c, s, &
-s, c /), &
(/2, 2/))
else if(what == ANALYTIC_EIGENVALUES) then
mat = reshape((/ 1.0_rk8, 0.0_rk8, &
0.0_rk8, am /), &
(/2, 2/))
else
!assert(0)
end if
! write(*,*) "calc value, elem: ", element, ", mat: ", mod(ii,2), mod(jj,2), mat(mod(ii,2), mod(jj,2)), "am ", am
! write(*,*) " matrix mat", mat
element = element * mat(mod(ii,2) + 1, mod(jj,2) + 1)
ii = ii / 2
jj = jj / 2
end do
!write(*,*) "returning value ", element
end function
function get_a(n) result (a)
implicit none
integer(kind=ik), intent(in) :: n
real(kind=rk8) :: a
real(kind=rk8), parameter :: largest_ev = 2.0
a = exp(log(largest_ev)/(2**n-1))
end function
end module
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment