Commit 214ddc43 authored by Andreas Marek's avatar Andreas Marek
Browse files

Fix single precision

parent 0d00252b
......@@ -282,7 +282,7 @@ print(" " + " \\\n ".join([
print("endif")
name = "test_skewsymmetric_real_double"
print("check_SCRIPTS += " + name)
print("check_SCRIPTS += " + name + "_extended.sh")
print("noinst_PROGRAMS += " + name)
print(name + "_SOURCES = test/Fortran/test_skewsymmetric.F90")
print(name + "_LDADD = $(test_program_ldadd)")
......@@ -293,7 +293,7 @@ print(" " + " \\\n ".join([
name = "test_skewsymmetric_real_single"
print("if WANT_SINGLE_PRECISION_REAL")
print("check_SCRIPTS += " + name)
print("check_SCRIPTS += " + name + "_extended.sh")
print("noinst_PROGRAMS += " + name)
print(name + "_SOURCES = test/Fortran/test_skewsymmetric.F90")
print(name + "_LDADD = $(test_program_ldadd)")
......
......@@ -175,7 +175,7 @@
#define PRECISION_NRM2 SNRM2
#define PRECISION_LASET SLASET
#define PRECISION_GER SGER
#define PRECISION_SCALSSSCAL
#define PRECISION_SCAL SSCAL
#define PRECISION_COPY SCOPY
#define PRECISION_AXPY SAXPY
#define cublas_PRECISION_GEMM cublas_SGEMM
......
......@@ -58,12 +58,16 @@ error: define exactly one of TEST_SINGLE or TEST_DOUBLE
#ifdef TEST_SINGLE
# define EV_TYPE real(kind=C_FLOAT)
# define EV_TYPE_COMPLEX complex(kind=C_FLOAT_COMPLEX)
# define MATRIX_TYPE_COMPLEX complex(kind=C_FLOAT_COMPLEX)
# ifdef TEST_REAL
# define MATRIX_TYPE real(kind=C_FLOAT)
# else
# define MATRIX_TYPE complex(kind=C_FLOAT_COMPLEX)
# endif
#else
# define MATRIX_TYPE_COMPLEX complex(kind=C_DOUBLE_COMPLEX)
# define EV_TYPE_COMPLEX complex(kind=C_DOUBLE_COMPLEX)
# define EV_TYPE real(kind=C_DOUBLE)
# ifdef TEST_REAL
# define MATRIX_TYPE real(kind=C_DOUBLE)
......@@ -72,14 +76,6 @@ error: define exactly one of TEST_SINGLE or TEST_DOUBLE
# endif
#endif
#ifdef TEST_SINGLE
#define MATRIX_TYPE_COMPLEX complex(kind=C_FLOAT_COMPLEX)
#define EV_TYPE_COMPLEX complex(kind=C_FLOAT_COMPLEX)
#else
#define MATRIX_TYPE_COMPLEX complex(kind=C_DOUBLE_COMPLEX)
#define EV_TYPE_COMPLEX complex(kind=C_DOUBLE_COMPLEX)
#endif
#ifdef TEST_REAL
# define AUTOTUNE_DOMAIN ELPA_AUTOTUNE_DOMAIN_REAL
#else
......@@ -199,11 +195,18 @@ program test
do j=1, na_cols
do i=1,na_rows
a_complex(i,j) = dcmplx(0.0, a_skewsymmetric(i,j))
enddo
do i=1,na_rows
#ifdef TEST_DOUBLE
a_complex(i,j) = dcmplx(0.0, a_skewsymmetric(i,j))
#endif
#ifdef TEST_SINGLE
a_complex(i,j) = cmplx(0.0, a_skewsymmetric(i,j))
#endif
enddo
enddo
z_complex(1:na_rows,1:na_cols) = a_complex(1:na_rows,1:na_cols)
as_complex(1:na_rows,1:na_cols) = a_complex(1:na_rows,1:na_cols)
......@@ -239,8 +242,8 @@ program test
status = check_correctness_evp_numeric_residuals_complex_double(na, nev, as_complex, z_complex, ev_complex, sc_desc, &
nblk, myid, np_rows,np_cols, my_prow, my_pcol)
#endif
! status = 0
! call check_status(status, myid)
status = 0
call check_status(status, myid)
#ifdef WITH_MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
......@@ -271,27 +274,36 @@ program test
call e_skewsymmetric%print_times("eigenvectors: skewsymmetric")
endif
! check eigenvalues
do i=1, na
if (myid == 0) then
! print *,"ev(", i,")=",ev_skewsymmetric(i)
#ifdef TEST_DOUBLE
if (abs(ev_complex(i)-ev_skewsymmetric(i))/abs(ev_complex(i)) .gt. 1e-10) then
#endif
#ifdef TEST_SINGLE
if (abs(ev_complex(i)-ev_skewsymmetric(i))/abs(ev_complex(i)) .gt. 1e-4) then
#endif
print *,"ev: i=",i,ev_complex(i),ev_skewsymmetric(i)
status = 1
endif
endif
enddo
! call check_status(status, myid)
z_complex(:,:) = 0
do j=1, na_cols
do i=1,na_rows
#ifdef TEST_DOUBLE
z_complex(i,j) = dcmplx(z_skewsymmetric(i,j), z_skewsymmetric(i,na_cols+j))
#endif
#ifdef TEST_SINGLE
z_complex(i,j) = cmplx(z_skewsymmetric(i,j), z_skewsymmetric(i,na_cols+j))
#endif
enddo
enddo
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
status = check_correctness_evp_numeric_residuals_ss(na, nev, as_skewsymmetric, z_complex, ev_skewsymmetric, &
sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
......
......@@ -60,15 +60,6 @@
complex(kind=rck), allocatable :: as_complex(:,:)
#ifndef WITH_MPI
#if REALCASE == 1
real(kind=rck) :: dnrm2, snrm2
#endif
#if COMPLEXCASE == 1
complex(kind=rck) :: zdotc, cdotc
#endif /* COMPLEXCASE */
#endif
integer(kind=ik) :: sc_desc(:)
integer(kind=ik) :: i, j, rowLocal, colLocal
......@@ -104,7 +95,6 @@
status = 0
! Setup complex matrices and eigenvalues
na_rows = size(as,dim=1)
na_cols = size(as,dim=2)
......@@ -131,9 +121,17 @@
xc = cmplx(0.0_rk,ev(i))
#endif
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call pzscal(na, xc, tmp1, 1, i, sc_desc, 1)
#else
call pcscal(na, xc, tmp1, 1, i, sc_desc, 1)
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
call zscal(na,xc,tmp1(:,i),1)
#else
call cscal(na,xc,tmp1(:,i),1)
#endif
#endif /* WITH_MPI */
enddo
......@@ -143,10 +141,19 @@
! tmp1 = A * Z
! as is original stored matrix, Z are the EVs
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call PZGEMM('N', 'N', na, nev, na, CONE, as_complex, 1, 1, sc_desc, &
z, 1, 1, sc_desc, CZERO, tmp1, 1, 1, sc_desc)
#else
call PCGEMM('N', 'N', na, nev, na, CONE, as_complex, 1, 1, sc_desc, &
z, 1, 1, sc_desc, CZERO, tmp1, 1, 1, sc_desc)
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
call ZGEMM('N','N',na,nev,na,CONE,as_complex,na,z,na,CZERO,tmp1,na)
#else
call CGEMM('N','N',na,nev,na,CONE,as_complex,na,z,na,CZERO,tmp1,na)
#endif
#endif /* WITH_MPI */
! tmp1 = A*Zi - Zi*EVi
......@@ -158,9 +165,17 @@
do i=1,nev
xc = (0.0_rk,0.0_rk)
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call PZDOTC(na, xc, tmp1, 1, i, sc_desc, 1, tmp1, 1, i, sc_desc, 1)
#else
call PCDOTC(na, xc, tmp1, 1, i, sc_desc, 1, tmp1, 1, i, sc_desc, 1)
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
xc = ZDOTC(na,tmp1,1,tmp1,1)
#else
xc = CDOTC(na,tmp1,1,tmp1,1)
#endif
#endif /* WITH_MPI */
errmax = max(errmax, sqrt(real(xc,kind=REAL_DATATYPE)))
enddo
......@@ -188,10 +203,20 @@
tmp2(:,:) = z(:,:)
tmp1 = 0
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call PZGEMM('C', 'N', nev, nev, na, CONE, z, 1, 1, &
sc_desc, tmp2, 1, 1, sc_desc, CZERO, tmp1, 1, 1, sc_desc)
#else
call PCGEMM('C', 'N', nev, nev, na, CONE, z, 1, 1, &
sc_desc, tmp2, 1, 1, sc_desc, CZERO, tmp1, 1, 1, sc_desc)
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
call ZGEMM('C','N',nev,nev,na,CONE,z,na,tmp2,na,CZERO,tmp1,na)
#else
call CGEMM('C','N',nev,nev,na,CONE,z,na,tmp2,na,CZERO,tmp1,na)
#endif
#endif /* WITH_MPI */
! First check, whether the elements on diagonal are 1 .. "normality" of the vectors
err = 0.0_rk
......@@ -211,9 +236,17 @@
! Initialize tmp2 to unit matrix
tmp2 = 0
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call PZLASET('A', nev, nev, CZERO, CONE, tmp2, 1, 1, sc_desc)
#else
call PCLASET('A', nev, nev, CZERO, CONE, tmp2, 1, 1, sc_desc)
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
call ZLASET('A',nev,nev,CZERO,CONE,tmp2,na)
#else
call CLASET('A',nev,nev,CZERO,CONE,tmp2,na)
#endif
#endif /* WITH_MPI */
! ! tmp1 = Z**T * Z - Unit Matrix
......@@ -286,12 +319,6 @@ function check_correctness_evp_numeric_residuals_&
MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
MATH_DATATYPE(kind=rck) :: xc
#ifndef WITH_MPI
#if COMPLEXCASE == 1
complex(kind=rck) :: zdotc, cdotc
#endif /* COMPLEXCASE */
#endif
integer(kind=ik) :: sc_desc(:)
integer(kind=ik) :: i, rowLocal, colLocal
......
Supports Markdown
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