Skip to content

GitLab

elpa
elpa
Commit 15a6cbc3 authored by Pavel Kus's avatar Pavel Kus Committed by Andreas Marek
Browse files
Options

templating test_analytic 

test analytic has been transformed to template to allow all
real/complex and single/double variants. However, at this
commit, only real double and real single variants are enabled
parent 7e701a7c
Hide whitespace changes
Inline Side-by-side
Showing with 395 additions and 223 deletions
Makefile.am
View file @ 15a6cbc3
......@@ -642,9 +642,8 @@ EXTRA_DIST = \
test/Fortran/assert.h \
test/Fortran/elpa_print_headers.F90 \
test/shared/test_check_correctness_template.F90 \
test/shared/test_check_correctness_template.F90 \
test/shared/test_prepare_matrix_template.F90 \
test/shared/test_prepare_matrix_template.F90 \
test/shared/test_analytic_template.F90 \
test_project/Makefile.am \
test_project/autogen.sh \
test_project/configure.ac \
......
generate_automake_test_programs.py
View file @ 15a6cbc3
......@@ -47,8 +47,7 @@ for m, g, t, p, d, s, l in product(
sorted(solver_flag.keys()),
sorted(layout_flag.keys())):
#todo: decide what tests we actually want
if(m == "analytic" and (g == 1 or t != "eigenvectors" or p == "single" or d == "complex")):
if(m == "analytic" and (g == 1 or t != "eigenvectors" or d == "complex" )):
continue
if(s == "scalapack_all" and (g == 1 or t != "eigenvectors" or p == "single" or d == "complex" or m != "analytic")):
......
test/shared/test_analytic.F90
View file @ 15a6cbc3
......@@ -47,112 +47,36 @@
module test_analytic
use test_util
interface prepare_matrix_analytic
module procedure prepare_matrix_analytic_complex_double
module procedure prepare_matrix_analytic_real_double
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure prepare_matrix_analytic_real_single
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure prepare_matrix_analytic_complex_single
#endif
end interface
interface check_correctness_analytic
module procedure check_correctness_analytic_complex_double
module procedure check_correctness_analytic_real_double
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure check_correctness_analytic_real_single
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure check_correctness_analytic_complex_single
#endif
end interface
integer(kind=ik), parameter, private :: num_primes = 3
integer(kind=ik), parameter, private :: primes(num_primes) = (/2,3,5/)
real(kind=rk8), parameter, private :: ZERO = 0.0_rk8, ONE = 1.0_rk8
contains
#include "../../src/general/prow_pcol.F90"
#include "../../src/general/map_global_to_local.F90"
subroutine prepare_matrix_analytic_real_double (na, a, nblk, myid, np_rows, &
np_cols, my_prow, my_pcol)
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(num_primes)
! for debug only, do it systematicaly somehow ... unit tests
! call check_module_sanity(myid)
if(.not. decompose(na, levels)) then
if(myid == 0) then
print *, "Analytic test can be run only with matrix sizes of the form 2^n * 3^m * 5^o"
stop 1
end if
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(na, globI, globJ)
end if
end do
end do
integer(kind=ik), parameter, private :: ANALYTIC_MATRIX = 0
integer(kind=ik), parameter, private :: ANALYTIC_EIGENVECTORS = 1
integer(kind=ik), parameter, private :: ANALYTIC_EIGENVALUES = 2
end subroutine
function check_correctness_analytic_real_double (na, nev, ev, z, nblk, myid, np_rows, &
np_cols, my_prow, my_pcol) result(status)
implicit none
integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol
integer(kind=ik) :: status, mpierr
real(kind=rk8), intent(inout) :: z(:,:)
real(kind=rk8), intent(inout) :: ev(:)
integer(kind=ik) :: globI, globJ, locI, locJ, levels(num_primes)
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 = ZERO
max_ev_diff = ZERO
do globJ = 1, na
diff = abs(ev(globJ) - analytic_eigenvalues(na, globJ))
max_ev_diff = max(diff, max_ev_diff)
end do
do globJ = 1, nev
! calculated eigenvector can be in opposite direction
max_curr_z_diff_minus = ZERO
max_curr_z_diff_plus = ZERO
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(na, 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 *, 'Maximum error in eigenvalues :', max_ev_diff
if(myid == 0) print *, 'Maximum error in eigenvectors :', glob_max_z_diff
status = 0
if (nev .gt. 2) then
if (max_ev_diff .gt. 5e-14_rk8 .or. max_ev_diff .eq. ZERO) status = 1
if (glob_max_z_diff .gt. 6e-11_rk8 .or. glob_max_z_diff .eq. ZERO) status = 1
else
if (max_ev_diff .gt. 5e-14_rk8) status = 1
if (glob_max_z_diff .gt. 6e-11_rk8) status = 1
endif
end function
contains
function decompose(num, decomposition) result(possible)
implicit none
......@@ -187,124 +111,6 @@ module test_analytic
end do
end function
#define ANALYTIC_MATRIX 0
#define ANALYTIC_EIGENVECTORS 1
#define ANALYTIC_EIGENVALUES 2
function analytic_matrix(na, i, j) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i, j
real(kind=rk8) :: element
element = analytic(na, i, j, ANALYTIC_MATRIX)
end function
function analytic_eigenvectors(na, i, j) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i, j
real(kind=rk8) :: element
element = analytic(na, i, j, ANALYTIC_EIGENVECTORS)
end function
function analytic_eigenvalues(na, i) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i
real(kind=rk8) :: element
element = analytic(na, i, i, ANALYTIC_EIGENVALUES)
end function
function analytic(na, i, j, what) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i, j, what
real(kind=rk8) :: element, am
real(kind=rk8) :: a, s, c, mat2x2(2,2), mat(5,5)
integer(kind=ik) :: levels(num_primes)
integer(kind=ik) :: ii, jj, m, prime_id, prime, total_level, level
real(kind=rk8), parameter :: largest_ev = 2.0_rk8
assert(i <= na)
assert(j <= na)
assert(i >= 0)
assert(j >= 0)
assert(decompose(na, levels))
! go to zero-based indexing
ii = i - 1
jj = j - 1
if (na .gt. 2) then
a = exp(log(largest_ev)/(na-1))
else
a = exp(log(largest_ev)/(1))
endif
s = 0.5_rk8
c = sqrt(3.0_rk8)/2.0_rk8
element = ONE
total_level = 0
am = a
do prime_id = 1,num_primes
prime = primes(prime_id)
do level = 1, levels(prime_id)
total_level = total_level + 1
if(what == ANALYTIC_MATRIX) then
mat2x2 = reshape((/ c*c + am**(prime-1) * s*s, (ONE-am**(prime-1)) * s*c, &
(ONE-am**(prime-1)) * s*c, s*s + am**(prime-1) * c*c /), &
(/2, 2/))
else if(what == ANALYTIC_EIGENVECTORS) then
mat2x2 = reshape((/ c, s, &
-s, c /), &
(/2, 2/))
else if(what == ANALYTIC_EIGENVALUES) then
mat2x2 = reshape((/ ONE, ZERO, &
ZERO, am**(prime-1) /), &
(/2, 2/))
else
assert(.false.)
end if
mat = ZERO
if(prime == 2) then
mat(1:2, 1:2) = mat2x2
else if(prime == 3) then
mat((/1,3/),(/1,3/)) = mat2x2
if(what == ANALYTIC_EIGENVECTORS) then
mat(2,2) = ONE
else
mat(2,2) = am
end if
else if(prime == 5) then
mat((/1,5/),(/1,5/)) = mat2x2
if(what == ANALYTIC_EIGENVECTORS) then
mat(2,2) = ONE
mat(3,3) = ONE
mat(4,4) = ONE
else
mat(2,2) = am
mat(3,3) = am**2
mat(4,4) = am**3
end if
else
assert(.false.)
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,prime) + 1, mod(jj,prime) + 1)
ii = ii / prime
jj = jj / prime
am = am**prime
end do
end do
!write(*,*) "returning value ", element
end function
subroutine print_matrix(myid, na, mat, mat_name)
implicit none
integer(kind=ik), intent(in) :: myid, na
......@@ -323,25 +129,70 @@ module test_analytic
write(*,*)
end subroutine
#include "../../src/general/prow_pcol.F90"
#include "../../src/general/map_global_to_local.F90"
#define COMPLEXCASE 1
#define DOUBLE_PRECISION 1
#include "../../src/general/precision_macros.h"
#include "test_analytic_template.F90"
#undef DOUBLE_PRECISION
#undef COMPLEXCASE
#ifdef WANT_SINGLE_PRECISION_COMPLEX
#define COMPLEXCASE 1
#define SINGLE_PRECISION 1
#include "../../src/general/precision_macros.h"
#include "test_analytic_template.F90"
#undef SINGLE_PRECISION
#undef COMPLEXCASE
#endif /* WANT_SINGLE_PRECISION_COMPLEX */
#define REALCASE 1
#define DOUBLE_PRECISION 1
#include "../../src/general/precision_macros.h"
#include "test_analytic_template.F90"
#undef DOUBLE_PRECISION
#undef REALCASE
#ifdef WANT_SINGLE_PRECISION_REAL
#define REALCASE 1
#define SINGLE_PRECISION 1
#include "../../src/general/precision_macros.h"
#include "test_analytic_template.F90"
#undef SINGLE_PRECISION
#undef REALCASE
#endif /* WANT_SINGLE_PRECISION_REAL */
subroutine check_matrices(myid, na)
implicit none
integer(kind=ik), intent(in) :: myid, na
real(kind=rk8) :: A(na, na), S(na, na), L(na, na)
real(kind=rk8) :: A(na, na), S(na, na), L(na, na), res(na, na)
integer(kind=ik) :: i, j, decomposition(num_primes)
assert(decompose(na, decomposition))
do i = 1, na
do j = 1, na
A(i,j) = analytic_matrix(na, i, j)
S(i,j) = analytic_eigenvectors(na, i, j)
L(i,j) = analytic(na, i, j, ANALYTIC_EIGENVALUES)
A(i,j) = analytic_matrix_real_double(na, i, j)
S(i,j) = analytic_eigenvectors_real_double(na, i, j)
L(i,j) = analytic_real_double(na, i, j, ANALYTIC_EIGENVALUES)
end do
end do
call print_matrix(myid, na, A, "A")
call print_matrix(myid, na, S, "S")
call print_matrix(myid, na, L, "L")
res = matmul(A,S) - matmul(S,L)
assert(maxval(abs(res)) < 1e-8)
!call print_matrix(myid, na, A, "A")
!call print_matrix(myid, na, S, "S")
!call print_matrix(myid, na, L, "L")
!call print_matrix(myid, na, res , "res")
end subroutine
......@@ -361,6 +212,8 @@ module test_analytic
call check_matrices(myid, 5)
call check_matrices(myid, 6)
call check_matrices(myid, 10)
call check_matrices(myid, 25)
call check_matrices(myid, 150)
if(myid == 0) print *, "Checking test_analytic module sanity.... DONE"
......
test/shared/test_analytic_template.F90 0 → 100644
View file @ 15a6cbc3
! (c) Copyright Pavel Kus, 2017, MPCDF
!
! This file is part of ELPA.
!
! The ELPA library was originally created by the ELPA consortium,
! consisting of the following organizations:
!
! - Max Planck Computing and Data Facility (MPCDF), formerly known as
! Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
! - Bergische Universität Wuppertal, Lehrstuhl für angewandte
! Informatik,
! - Technische Universität München, Lehrstuhl für Informatik mit
! Schwerpunkt Wissenschaftliches Rechnen ,
! - Fritz-Haber-Institut, Berlin, Abt. Theorie,
! - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
! Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
! and
! - IBM Deutschland GmbH
!
!
! More information can be found here:
! http://elpa.mpcdf.mpg.de/
!
! ELPA is free software: you can redistribute it and/or modify
! it under the terms of the version 3 of the license of the
! GNU Lesser General Public License as published by the Free
! Software Foundation.
!
! ELPA is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with ELPA. If not, see <http://www.gnu.org/licenses/>
!
! ELPA reflects a substantial effort on the part of the original
! ELPA consortium, and we ask you to respect the spirit of the
! license that we chose: i.e., please contribute any changes you
! may have back to the original ELPA library distribution, and keep
! any derivatives of ELPA under the same license that we chose for
! the original distribution, the GNU Lesser General Public License.
subroutine prepare_matrix_analytic_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
implicit none
integer(kind=ik), intent(in) :: na, nblk, myid, np_rows, np_cols, my_prow, my_pcol
MATH_DATATYPE(kind=REAL_DATATYPE), intent(inout) :: a(:,:)
integer(kind=ik) :: globI, globJ, locI, locJ, levels(num_primes)
! for debug only, do it systematicaly somehow ... unit tests
call check_module_sanity(myid)
if(.not. decompose(na, levels)) then
if(myid == 0) then
print *, "Analytic test can be run only with matrix sizes of the form 2^n * 3^m * 5^o"
stop 1
end if
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_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, globI, globJ)
end if
end do
end do
end subroutine
function check_correctness_analytic_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol) 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(:)
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_minus, max_curr_z_diff_plus
#ifdef DOUBLE_PRECISION
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
#endif
real(kind=rk) :: computed_ev, expected_ev
MATH_DATATYPE(kind=rck) :: computed_z, expected_z
if(.not. decompose(na, levels)) then
print *, "can not decomopse matrix size"
stop 1
end if
!call print_matrix(myid, na, z, "z")
max_z_diff = 0.0_rk
max_ev_diff = 0.0_rk
do globJ = 1, na
computed_ev = ev(globJ)
expected_ev = analytic_eigenvalues_real_&
&PRECISION&
&(na, globJ)
diff = abs(computed_ev - expected_ev)
max_ev_diff = max(diff, max_ev_diff)
end do
do globJ = 1, nev
! calculated eigenvector can be in opposite direction
max_curr_z_diff_minus = 0.0_rk
max_curr_z_diff_plus = 0.0_rk
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 = z(locI, locJ)
expected_z = analytic_eigenvectors_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, globI, globJ)
max_curr_z_diff_minus = max(abs(computed_z - expected_z), max_curr_z_diff_minus)
max_curr_z_diff_plus = max(abs(computed_z + expected_z), 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_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else
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
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
else
if (max_ev_diff .gt. tol_eigenvalues) status = 1
if (glob_max_z_diff .gt. tol_eigenvectors) status = 1
endif
end function
function analytic_matrix_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, i, j) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i, j
MATH_DATATYPE(kind=REAL_DATATYPE) :: element
element = analytic_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, i, j, ANALYTIC_MATRIX)
end function
function analytic_eigenvectors_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, i, j) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i, j
MATH_DATATYPE(kind=REAL_DATATYPE) :: element
element = analytic_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, i, j, ANALYTIC_EIGENVECTORS)
end function
function analytic_eigenvalues_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, i) result(element)
implicit none
integer(kind=ik), intent(in) :: na, i
real(kind=REAL_DATATYPE) :: element
element = analytic_real_&
&PRECISION&
&(na, i, i, ANALYTIC_EIGENVALUES)
end function
function analytic_&
&MATH_DATATYPE&
&_&
&PRECISION&
&(na, i, j, what) result(element)
implicit none
#include "../../src/general/precision_kinds.F90"
integer(kind=ik), intent(in) :: na, i, j, what
MATH_DATATYPE(kind=rck) :: element, mat2x2(2,2), mat(5,5)
real(kind=rk) :: a, am, amp
integer(kind=ik) :: levels(num_primes)
integer(kind=ik) :: ii, jj, m, prime_id, prime, total_level, level
real(kind=rk), parameter :: s = 0.5_rk
real(kind=rk), parameter :: c = 0.86602540378443864679_rk
real(kind=rk), parameter :: sq2 = 1.4142135623730950488_rk
real(kind=rk), parameter :: largest_ev = 2.0_rk
assert(i <= na)
assert(j <= na)
assert(i >= 0)
assert(j >= 0)
assert(decompose(na, levels))
! go to zero-based indexing
ii = i - 1
jj = j - 1
if (na .gt. 2) then
a = exp(log(largest_ev)/(na-1))
else