!#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