Commit 8237833c authored by Andreas Marek's avatar Andreas Marek
Browse files

Cleanup of check_correctness

parent 3ac5664b
......@@ -135,13 +135,7 @@ program test
EV_TYPE, allocatable :: ev(:)
#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
EV_TYPE, allocatable :: d(:), sd(:), ev_analytic(:), ds(:), sds(:)
EV_TYPE :: diagonalELement, subdiagonalElement, tmp, maxerr
#ifdef TEST_DOUBLE
EV_TYPE, parameter :: pi = 3.141592653589793238462643383279_rk8
#else
EV_TYPE, parameter :: pi = 3.1415926535897932_rk4
#endif
integer :: loctmp ,rowLocal, colLocal, j,ii
EV_TYPE :: diagonalELement, subdiagonalElement
#endif
......@@ -221,7 +215,17 @@ program test
#endif
#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
call prepare_toeplitz_matrix(na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
#ifdef TEST_SINGLE
diagonalElement = 0.45_c_float
subdiagonalElement = 0.78_c_float
#else
diagonalElement = 0.45_c_double
subdiagonalElement = 0.78_c_double
#endif
call prepare_toeplitz_matrix(na, diagonalElement, subdiagonalElement, &
d, sd, ds, sds, a, as, nblk, np_rows, &
np_cols, my_prow, my_pcol)
#endif
if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
......@@ -355,55 +359,9 @@ program test
if (myid == 0) print *, ""
#endif
#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
status = 0
! analytic solution
diagonalElement = ds(1)
subdiagonalElement = sds(1)
do ii=1, na
#ifdef TEST_DOUBLE
ev_analytic(ii) = diagonalElement + 2.0 * subdiagonalElement *cos( pi*real(ii,kind=rk8)/ real(na+1,kind=rk8) )
#else
ev_analytic(ii) = diagonalElement + 2.0 * subdiagonalElement *cos( pi*real(ii,kind=rk4)/ real(na+1,kind=rk4) )
#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 ii=2, na
tmp = ev_analytic(ii)
do j= ii, na
if (ev_analytic(j) .lt. tmp) then
tmp = ev_analytic(j)
loctmp = j
endif
enddo
ev_analytic(loctmp) = ev_analytic(ii)
ev_analytic(ii) = tmp
enddo
! compute a simple error max of eigenvalues
maxerr = 0.0
maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)
#ifdef TEST_DOUBLE
if (maxerr .gt. 8.e-13) then
#else
if (maxerr .gt. 8.e-4) then
#endif
status = 1
if (myid .eq. 0) then
print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
endif
endif
status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, &
subdiagonalElement, ev, z, myid)
if (status /= 0) then
call exit(status)
......
......@@ -56,6 +56,16 @@ module test_check_correctness
#endif
end interface
interface check_correctness_eigenvalues_toeplitz
module procedure check_correctness_eigenvalues_toeplitz_complex_double
module procedure check_correctness_eigenvalues_toeplitz_real_double
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure check_correctness_eigenvalues_toeplitz_real_single
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure check_correctness_eigenvalues_toeplitz_complex_single
#endif
end interface
contains
#define COMPLEXCASE 1
......
......@@ -433,5 +433,83 @@ function check_correctness_&
end function
function check_correctness_eigenvalues_toeplitz_&
&MATH_DATATYPE&
&_&
&PRECISION&
& (na, diagonalElement, subdiagonalElement, ev, z, myid) result(status)
use iso_c_binding
implicit none
integer :: status, ii, j, myid
integer, intent(in) :: na
real(kind=C_DATATYPE_KIND) :: diagonalElement, subdiagonalElement
real(kind=C_DATATYPE_KIND) :: ev_analytic(na), ev(na)
#if REALCASE == 1
real(kind=C_DATATYPE_KIND) :: z(:,:)
#endif
#if COMPLEXCASE == 1
complex(kind=C_DATATYPE_KIND) :: z(:,:)
#endif
#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
real(kind=C_DATATYPE_KIND), parameter :: pi = 3.141592653589793238462643383279_c_double
#else
real(kind=C_DATATYPE_KIND), parameter :: pi = 3.1415926535897932_c_float
#endif
real(kind=C_DATATYPE_KIND) :: tmp, maxerr
integer :: loctmp
status = 0
! analytic solution
do ii=1, na
#if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
ev_analytic(ii) = diagonalElement + 2.0_c_double * &
subdiagonalElement *cos( pi*real(ii,kind=c_double)/ &
real(na+1,kind=c_double) )
#else
ev_analytic(ii) = diagonalElement + 2.0_c_float * &
subdiagonalElement *cos( pi*real(ii,kind=c_float)/ &
real(na+1,kind=c_float) )
#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 ii=2, na
tmp = ev_analytic(ii)
do j= ii, na
if (ev_analytic(j) .lt. tmp) then
tmp = ev_analytic(j)
loctmp = j
endif
enddo
ev_analytic(loctmp) = ev_analytic(ii)
ev_analytic(ii) = 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 *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
endif
endif
end function
! vim: syntax=fortran
......@@ -210,8 +210,8 @@ subroutine prepare_matrix_&
&MATH_DATATYPE&
&_&
&PRECISION&
& (na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
use iso_c_binding
& (na, diagonalElement, subdiagonalElement, d, sd, ds, sds, a, as, &
nblk, np_rows, np_cols, my_prow, my_pcol)
use test_util
implicit none
......@@ -220,36 +220,14 @@ subroutine prepare_matrix_&
real(kind=C_DATATYPE_KIND) :: d(:), sd(:), ds(:), sds(:)
#if REALCASE == 1
real(kind=C_DATATYPE_KIND) :: z(:,:), a(:,:), as(:,:)
real(kind=C_DATATYPE_KIND) :: a(:,:), as(:,:)
#endif
#if COMPLEXCASE == 1
complex(kind=C_DATATYPE_KIND) :: z(:,:), a(:,:), as(:,:)
complex(kind=C_DATATYPE_KIND) :: a(:,:), as(:,:)
#endif
integer :: ii, rowLocal, colLocal
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
diagonalElement = 0.45_rk8
subdiagonalElement = 0.78_rk8
#else
diagonalElement = 0.45_rk4
subdiagonalElement = 0.78_rk4
#endif
#endif /* REALCASE */
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
diagonalElement = 0.45_rk8
subdiagonalElement = 0.78_rk8
#else
diagonalElement = 0.45_rk4
subdiagonalElement = 0.78_rk4
#endif
#endif /* COMPLEXCASE */
d(:) = diagonalElement
sd(:) = subdiagonalElement
......
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