Commit 44122b1f authored by Andreas Marek's avatar Andreas Marek

Frank matrix in test programs

parent 4c7521e0
......@@ -23,6 +23,13 @@ gpu_flag = {
matrix_flag = {
"random" : "-DTEST_MATRIX_RANDOM",
"analytic" : "-DTEST_MATRIX_ANALYTIC",
"toeplitz" : "-DTEST_MATRIX_TOEPLITZ",
"frank" : "-DTEST_MATRIX_FRANK",
}
qr_flag = {
0 : "-DTEST_QR_DECOMPOSITION=0",
1 : "-DTEST_QR_DECOMPOSITION=1",
}
test_type_flag = {
......@@ -31,7 +38,6 @@ test_type_flag = {
"solve_tridiagonal" : "-DTEST_SOLVE_TRIDIAGONAL",
"cholesky" : "-DTEST_CHOLESKY",
"hermitian_multiply" : "-DTEST_HERMITIAN_MULTIPLY",
"qr" : "-DTEST_QR_DECOMPOSITION",
}
layout_flag = {
......@@ -39,31 +45,49 @@ layout_flag = {
"square" : ""
}
for m, g, t, p, d, s, l in product(
for m, g, q, t, p, d, s, l in product(
sorted(matrix_flag.keys()),
sorted(gpu_flag.keys()),
sorted(qr_flag.keys()),
sorted(test_type_flag.keys()),
sorted(prec_flag.keys()),
sorted(domain_flag.keys()),
sorted(solver_flag.keys()),
sorted(layout_flag.keys())):
# exclude some test combinations
# analytic tests only for "eigenvectors" and not on GPU
if(m == "analytic" and (g == 1 or t != "eigenvectors")):
continue
# Frank tests only for "eigenvectors" and eigenvalues and real double precision case
if(m == "frank" and ((t != "eigenvectors" or t != "eigenvalues") and (d !="real" or p !="double"))):
continue
if(s in ["scalapack_all", "scalapack_part"] and (g == 1 or t != "eigenvectors" or m != "analytic")):
continue
if (t == "solve_tridiagonal" and (s == "2stage" or d == "complex")):
# solve tridiagonal only for real toeplitz matrix in 1stage
if (t == "solve_tridiagonal" and (s != "1stage" or d !="real" or m != "toeplitz")):
continue
if (t == "cholesky" and (s == "2stage")):
# cholesky tests only 1stage and teoplitz matrix
if (t == "cholesky" and (m != "toeplitz" or s == "2stage")):
continue
if (t == "eigenvalues" and (m == "random")):
continue
if (t == "hermitian_multiply" and (s == "2stage")):
continue
if (t == "qr" and (s == "1stage" or d == "complex")):
if (t == "hermitian_multiply" and (m == "toeplitz")):
continue
# qr only for 2stage real
if (q == 1 and (s != "2stage" or d != "real" or t != "eigenvectors" or g == 1 or m != "random")):
continue
for kernel in ["all_kernels", "default_kernel"] if s == "2stage" else ["nokernel"]:
......@@ -102,13 +126,18 @@ for m, g, t, p, d, s, l in product(
raise Exception("Oh no!")
endifs += 1
name = "test_{0}_{1}_{2}_{3}{4}{5}{6}{7}".format(
name = "test_{0}_{1}_{2}_{3}{4}_{5}{6}{7}{8}".format(
d, p, t, s,
"" if kernel == "nokernel" else "_" + kernel,
"_gpu" if g else "",
"_analytic" if m == "analytic" else "",
"gpu_" if g else "",
"qr_" if q else "",
m,
"_all_layouts" if l == "all_layouts" else "")
print("if BUILD_KCOMPUTER")
print("bin_PROGRAMS += " + name)
print("else")
print("noinst_PROGRAMS += " + name)
print("endif")
print("check_SCRIPTS += " + name + ".sh")
print(name + "_SOURCES = test/Fortran/test.F90")
print(name + "_LDADD = $(test_program_ldadd)")
......@@ -120,6 +149,7 @@ for m, g, t, p, d, s, l in product(
test_type_flag[t],
solver_flag[s],
gpu_flag[g],
qr_flag[q],
matrix_flag[m]] + extra_flags))
print("endif\n" * endifs)
This diff is collapsed.
......@@ -85,32 +85,39 @@
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals) result(status)
&(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals, check_eigenvectors) result(status)
implicit none
#include "../../src/general/precision_kinds.F90"
integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol
integer(kind=ik) :: status, mpierr
MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:)
real(kind=rk), intent(inout) :: ev(:)
logical, intent(in) :: check_all_evals
integer(kind=ik) :: globI, globJ, locI, locJ, levels(num_primes)
real(kind=rk) :: diff, max_z_diff, max_ev_diff, glob_max_z_diff, max_curr_z_diff
integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, &
np_cols, my_prow, my_pcol
integer(kind=ik) :: status, mpierr
MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:)
real(kind=rk), intent(inout) :: ev(:)
logical, intent(in) :: check_all_evals, check_eigenvectors
integer(kind=ik) :: globI, globJ, locI, locJ, &
levels(num_primes)
real(kind=rk) :: diff, max_z_diff, max_ev_diff, &
glob_max_z_diff, max_curr_z_diff
#ifdef DOUBLE_PRECISION
real(kind=rk), parameter :: tol_eigenvalues = 5e-14_rk8
real(kind=rk), parameter :: tol_eigenvectors = 6e-11_rk8
real(kind=rk), parameter :: tol_eigenvalues = 5e-14_rk8
real(kind=rk), parameter :: tol_eigenvectors = 6e-11_rk8
#endif
#ifdef SINGLE_PRECISION
! tolerance needs to be very high due to qr tests
! it should be distinguished somehow!
real(kind=rk), parameter :: tol_eigenvalues = 7e-6_rk4
real(kind=rk), parameter :: tol_eigenvectors = 4e-3_rk4
real(kind=rk), parameter :: tol_eigenvalues = 7e-6_rk4
real(kind=rk), parameter :: tol_eigenvectors = 4e-3_rk4
#endif
real(kind=rk) :: computed_ev, expected_ev
MATH_DATATYPE(kind=rck) :: computed_z, expected_z
real(kind=rk) :: computed_ev, expected_ev
MATH_DATATYPE(kind=rck) :: computed_z, expected_z
MATH_DATATYPE(kind=rck) :: max_value_for_normalization, computed_z_on_max_position, normalization_quotient
integer(kind=ik) :: max_value_idx, rank_with_max, rank_with_max_reduced, num_checked_evals
MATH_DATATYPE(kind=rck) :: max_value_for_normalization, &
computed_z_on_max_position, &
normalization_quotient
integer(kind=ik) :: max_value_idx, rank_with_max, &
rank_with_max_reduced, &
num_checked_evals
if(.not. decompose(na, levels)) then
......@@ -193,14 +200,21 @@
glob_max_z_diff = max_z_diff
#endif
if(myid == 0) print *, 'Maximum error in eigenvalues :', max_ev_diff
if(myid == 0) print *, 'Maximum error in eigenvectors :', glob_max_z_diff
if (check_eigenvectors) then
if(myid == 0) print *, 'Maximum error in eigenvectors :', glob_max_z_diff
endif
status = 0
if (nev .gt. 2) then
if (max_ev_diff .gt. tol_eigenvalues .or. max_ev_diff .eq. 0.0_rk) status = 1
if (glob_max_z_diff .gt. tol_eigenvectors .or. glob_max_z_diff .eq. 0.0_rk) status = 1
if (check_eigenvectors) then
if (glob_max_z_diff .gt. tol_eigenvectors .or. glob_max_z_diff .eq. 0.0_rk) status = 1
endif
else
if (max_ev_diff .gt. tol_eigenvalues) status = 1
if (glob_max_z_diff .gt. tol_eigenvectors) status = 1
if (check_eigenvectors) then
if (glob_max_z_diff .gt. tol_eigenvectors) status = 1
endif
endif
end function
......
......@@ -67,6 +67,17 @@ module test_check_correctness
#endif
end interface
interface check_correctness_eigenvalues_frank
module procedure check_correctness_eigenvalues_frank_complex_double
module procedure check_correctness_eigenvalues_frank_real_double
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure check_correctness_eigenvalues_frank_real_single
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure check_correctness_eigenvalues_frank_complex_single
#endif
end interface
interface check_correctness_cholesky
module procedure check_correctness_cholesky_complex_double
module procedure check_correctness_cholesky_real_double
......
......@@ -192,17 +192,17 @@
#if REALCASE == 1
if (nev .ge. 2) then
#ifdef DOUBLE_PRECISION_REAL
if (errmax .gt. 5e-12_rk8 .or. errmax .eq. 0.0_rk8) then
if (errmax .gt. 5e-11_rk8 .or. errmax .eq. 0.0_rk8) then
#else
if (errmax .gt. 3e-3_rk4 .or. errmax .eq. 0.0_rk4) then
if (errmax .gt. 3e-2_rk4 .or. errmax .eq. 0.0_rk4) then
#endif
status = 1
endif
else
#ifdef DOUBLE_PRECISION_REAL
if (errmax .gt. 5e-12_rk8) then
if (errmax .gt. 5e-11_rk8) then
#else
if (errmax .gt. 3e-3_rk4) then
if (errmax .gt. 2e-2_rk4) then
#endif
status = 1
endif
......@@ -211,17 +211,17 @@
#if COMPLEXCASE == 1
if (nev .gt. 2) then
#ifdef DOUBLE_PRECISION_COMPLEX
if (errmax .gt. 5e-12_rk8 .or. errmax .eq. 0.0_rk8) then
if (errmax .gt. 5e-11_rk8 .or. errmax .eq. 0.0_rk8) then
#else
if (errmax .gt. 3e-3_rk4 .or. errmax .eq. 0.0_rk4) then
if (errmax .gt. 3e-2_rk4 .or. errmax .eq. 0.0_rk4) then
#endif
status =1
endif
else
#ifdef DOUBLE_PRECISION_COMPLEX
if (errmax .gt. 5e-12_rk8) then
if (errmax .gt. 5e-11_rk8) then
#else
if (errmax .gt. 3e-3_rk4) then
if (errmax .gt. 3e-2_rk4) then
#endif
status = 1
endif
......@@ -289,9 +289,9 @@
endif
else
#ifdef DOUBLE_PRECISION_REAL
if (errmax .gt. 5e-12_rk8) then
if (errmax .gt. 5e-11_rk8) then
#else
if (errmax .gt. 9e-4_rk4) then
if (errmax .gt. 9e-2_rk4) then
#endif
status = 1
endif
......@@ -308,9 +308,9 @@
endif
else
#ifdef DOUBLE_PRECISION_COMPLEX
if (errmax .gt. 5e-12_rk8) then
if (errmax .gt. 5e-11_rk8) then
#else
if (errmax .gt. 9e-4_rk4) then
if (errmax .gt. 9e-3_rk4) then
#endif
status = 1
endif
......@@ -449,6 +449,7 @@ function check_correctness_evp_numeric_residuals_&
#endif
status = 1
if (myid .eq. 0) then
print *,"Result of Toeplitz matrix test: "
print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
endif
endif
......@@ -756,7 +757,7 @@ function check_correctness_evp_numeric_residuals_&
#ifdef DOUBLE_PRECISION_COMPLEX
complex(kind=ck8), parameter :: CZERO = (0.0_rk8,0.0_rk8), CONE = (1.0_rk8,0.0_rk8)
#else
complex(kind=ck4), parameter :: CZERO = (0.0_rk4,0.0_rk4), CONE = (1.0_rk4,0.0_rk8)
complex(kind=ck4), parameter :: CZERO = (0.0_rk4,0.0_rk4), CONE = (1.0_rk4,0.0_rk4)
#endif
#ifdef WITH_MPI
......@@ -966,7 +967,80 @@ function check_correctness_evp_numeric_residuals_&
#endif
end function
function check_correctness_eigenvalues_frank_&
&MATH_DATATYPE&
&_&
&PRECISION&
& (na, ev, z, myid) result(status)
use iso_c_binding
implicit none
#include "../../src/general/precision_kinds.F90"
integer :: status, i, j, myid
integer, intent(in) :: na
real(kind=rck) :: ev_analytic(na), ev(na)
MATH_DATATYPE(kind=rck) :: z(:,:)
#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
real(kind=rck), parameter :: pi = 3.141592653589793238462643383279_c_double
#else
real(kind=rck), parameter :: pi = 3.1415926535897932_c_float
#endif
real(kind=rck) :: tmp, maxerr
integer :: loctmp
status = 0
! analytic solution
do i = 1, na
j = na - i
#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
ev_analytic(i) = pi * (2.0_c_double * real(j,kind=c_double) + 1.0_c_double) / &
(2.0_c_double * real(na,kind=c_double) + 1.0_c_double)
ev_analytic(i) = 0.5_c_double / (1.0_c_double - cos(ev_analytic(i)))
#else
ev_analytic(i) = pi * (2.0_c_float * real(j,kind=c_float) + 1.0_c_float) / &
(2.0_c_float * real(na,kind=c_float) + 1.0_c_float)
ev_analytic(i) = 0.5_c_float / (1.0_c_float - cos(ev_analytic(i)))
#endif
enddo
! sort analytic solution:
! this hack is neither elegant, nor optimized: for huge matrixes it might be expensive
! a proper sorting algorithmus might be implemented here
tmp = minval(ev_analytic)
loctmp = minloc(ev_analytic, 1)
ev_analytic(loctmp) = ev_analytic(1)
ev_analytic(1) = tmp
do i=2, na
tmp = ev_analytic(i)
do j= i, na
if (ev_analytic(j) .lt. tmp) then
tmp = ev_analytic(j)
loctmp = j
endif
enddo
ev_analytic(loctmp) = ev_analytic(i)
ev_analytic(i) = tmp
enddo
! compute a simple error max of eigenvalues
maxerr = 0.0
maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)
#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
if (maxerr .gt. 8.e-13_c_double) then
#else
if (maxerr .gt. 8.e-4_c_float) then
#endif
status = 1
if (myid .eq. 0) then
print *,"Result of Frank matrix test: "
print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
endif
endif
end function
! vim: syntax=fortran
......@@ -69,6 +69,19 @@ module test_prepare_matrix
#endif
end interface
interface prepare_matrix_frank
module procedure prepare_matrix_frank_complex_double
module procedure prepare_matrix_frank_real_double
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure prepare_matrix_frank_real_single
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure prepare_matrix_frank_complex_single
#endif
end interface
private prows, pcols, map_global_array_index_to_local_index
contains
......
......@@ -320,6 +320,41 @@ subroutine prepare_matrix_random_&
#endif
end subroutine
subroutine prepare_matrix_frank_&
&MATH_DATATYPE&
&_&
&PRECISION&
& (na, a, z, as, nblk, np_rows, np_cols, my_prow, my_pcol)
use test_util
implicit none
integer, intent(in) :: na, nblk, np_rows, np_cols, my_prow, my_pcol
#if REALCASE == 1
real(kind=C_DATATYPE_KIND) :: a(:,:), z(:,:), as(:,:)
#endif
#if COMPLEXCASE == 1
complex(kind=C_DATATYPE_KIND) :: a(:,:), z(:,:), as(:,:)
#endif
integer :: i, j, rowLocal, colLocal
do i = 1, na
do j = 1, na
if (map_global_array_index_to_local_index(i, j, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
if (j .le. i) then
a(rowLocal,colLocal) = real((na+1-i), kind=C_DATATYPE_KIND) / real(na, kind=C_DATATYPE_KIND)
else
a(rowLocal,colLocal) = real((na+1-j), kind=C_DATATYPE_KIND) / real(na, kind=C_DATATYPE_KIND)
endif
endif
enddo
enddo
z(:,:) = a(:,:)
as(:,:) = a(:,:)
end subroutine
! vim: syntax=fortran
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