Commit bd683040 authored by Pavel Kus's avatar Pavel Kus
Browse files

unification in check_correctness_hermitian_multiply

parent d73f4721
......@@ -520,222 +520,78 @@ function check_correctness_evp_numeric_residuals_&
#include "../../src/general/precision_kinds.F90"
integer(kind=ik) :: status
integer(kind=ik), intent(in) :: na, myid, na_rows
#if REALCASE == 1
real(kind=rck), intent(in) :: a(:,:), b(:,:), c(:,:)
real(kind=rck), dimension(size(a,dim=1),size(a,dim=2)) :: tmp1, tmp2
MATH_DATATYPE(kind=rck), intent(in) :: a(:,:), b(:,:), c(:,:)
MATH_DATATYPE(kind=rck), dimension(size(a,dim=1),size(a,dim=2)) :: tmp1, tmp2
real(kind=rck) :: norm, normmax
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
real(kind=rck) :: pdlange
#else
real(kind=rck) :: pslange
#endif
real(kind=rck) :: p&
&BLAS_CHAR&
&lange
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
real(kind=rck) :: dlange
#else
real(kind=rck) :: slange
#endif
real(kind=rck) :: BLAS_CHAR&
&lange
#endif /* WITH_MPI */
#endif /* REALCASE */
#if COMPLEXCASE == 1
complex(kind=rck), intent(in) :: a(:,:), b(:,:), c(:,:)
complex(kind=rck), dimension(size(a,dim=1),size(a,dim=2)) :: tmp1, tmp2
real(kind=rck) :: norm, normmax
#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)
#endif
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_COMPLEX
real(kind=rck) :: pzlange
#else
real(kind=rck) :: pclange
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_COMPLEX
real(kind=rck) :: zlange
#else
real(kind=rck) :: clange
#endif
#endif /* WITH_MPI */
#endif /* COMPLEXCASE */
integer(kind=ik) :: sc_desc(:)
real(kind=rck) :: err, errmax
integer :: mpierr
status = 0
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
tmp1(:,:) = 0.0_rk8
#else
tmp1(:,:) = 0.0_rk4
#endif
#endif /* REALCASE */
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
tmp1(:,:) = (0.0_c_double, 0.0_c_double)
#else
tmp1(:,:) = (0.0_c_float, 0.0_c_float)
#endif
#endif /* COMPLEXCASE */
tmp1(:,:) = ZERO
#if REALCASE == 1
! tmp1 = a**T
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call pdtran(na, na, 1.0_rk8, a, 1, 1, sc_desc, 0.0_rk8, tmp1, 1, 1, sc_desc)
#else
call pstran(na, na, 1.0_rk4, a, 1, 1, sc_desc, 0.0_rk4, tmp1, 1, 1, sc_desc)
#endif
call p&
&BLAS_CHAR&
&tran(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
#else /* WITH_MPI */
tmp1 = transpose(a)
#endif /* WITH_MPI */
! tmp2 = tmp1 * b
#ifdef DOUBLE_PRECISION_REAL
#ifdef WITH_MPI
call pdgemm("N","N", na, na, na, 1.0_rk8, tmp1, 1, 1, sc_desc, b, 1, 1, &
sc_desc, 0.0_rk8, tmp2, 1, 1, sc_desc)
#else
call dgemm("N","N", na, na, na, 1.0_rk8, tmp1, na, b, na, 0.0_rk8, tmp2, na)
#endif
#else /* DOUBLE_PRECISION_REAL */
#ifdef WITH_MPI
call psgemm("N","N", na, na, na, 1.0_rk4, tmp1, 1, 1, sc_desc, b, 1, 1, &
sc_desc, 0.0_rk4, tmp2, 1, 1, sc_desc)
#else
call sgemm("N","N", na, na, na, 1.0_rk4, tmp1, na, b, na, 0.0_rk4, tmp2, na)
#endif
#endif /* DOUBLE_PRECISION_REAL */
#endif /* REALCASE == 1 */
#if COMPLEXCASE == 1
! tmp1 = a**H
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_COMPLEX
call pztranc(na, na, CONE, a, 1, 1, sc_desc, CZERO, tmp1, 1, 1, sc_desc)
#else
call pctranc(na, na, CONE, a, 1, 1, sc_desc, CZERO, tmp1, 1, 1, sc_desc)
#endif
call p&
&BLAS_CHAR&
&tranc(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
#else /* WITH_MPI */
tmp1 = transpose(conjg(a))
#endif /* WITH_MPI */
#endif /* COMPLEXCASE == 1 */
! tmp2 = tmp1 * b
#ifdef DOUBLE_PRECISION_COMPLEX
#ifdef WITH_MPI
call pzgemm("N","N", na, na, na, CONE, tmp1, 1, 1, sc_desc, b, 1, 1, &
sc_desc, CZERO, tmp2, 1, 1, sc_desc)
#else
call zgemm("N","N", na, na, na, CONE, tmp1, na, b, na, CZERO, tmp2, na)
#endif
#else /* DOUBLE_PRECISION_COMPLEX */
#ifdef WITH_MPI
call pcgemm("N","N", na, na, na, CONE, tmp1, 1, 1, sc_desc, b, 1, 1, &
sc_desc, CZERO, tmp2, 1, 1, sc_desc)
call p&
&BLAS_CHAR&
&gemm("N","N", na, na, na, ONE, tmp1, 1, 1, sc_desc, b, 1, 1, &
sc_desc, ZERO, tmp2, 1, 1, sc_desc)
#else
call cgemm("N","N", na, na, na, CONE, tmp1, na, b, na, CZERO, tmp2, na)
call BLAS_CHAR&
&gemm("N","N", na, na, na, ONE, tmp1, na, b, na, ZERO, tmp2, na)
#endif
#endif /* DOUBLE_PRECISION_COMPLEX */
#endif /* COMPLEXCASE == 1 */
! compare tmp2 with c
tmp2(:,:) = tmp2(:,:) - c(:,:)
#if REALCASE == 1
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
norm = pdlange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
#else
norm = pslange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
#endif
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_REAL
norm = dlange("M", na, na, tmp2, na_rows, tmp1)
#else
norm = slange("M", na, na, tmp2, na_rows, tmp1)
#endif
#endif /* WITH_MPI */
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_REAL
call mpi_allreduce(norm,normmax,1,MPI_REAL8,MPI_MAX,MPI_COMM_WORLD,mpierr)
#else
call mpi_allreduce(norm,normmax,1,MPI_REAL4,MPI_MAX,MPI_COMM_WORLD,mpierr)
#endif
#else /* WITH_MPI */
normmax = norm
#endif /* WITH_MPI */
#endif /* REALCASE == 1 */
#if COMPLEXCASE == 1
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_COMPLEX
norm = pzlange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
#else
norm = pclange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
#endif
norm = p&
&BLAS_CHAR&
&lange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
#else /* WITH_MPI */
#ifdef DOUBLE_PRECISION_COMPLEX
norm = zlange("M", na, na, tmp2, na_rows, tmp1)
#else
norm = clange("M", na, na, tmp2, na_rows, tmp1)
#endif
norm = BLAS_CHAR&
&lange("M", na, na, tmp2, na_rows, tmp1)
#endif /* WITH_MPI */
#ifdef WITH_MPI
#ifdef DOUBLE_PRECISION_COMPLEX
call mpi_allreduce(norm,normmax,1,MPI_REAL8,MPI_MAX,MPI_COMM_WORLD,mpierr)
#else
call mpi_allreduce(norm,normmax,1,MPI_REAL4,MPI_MAX,MPI_COMM_WORLD,mpierr)
#endif
call mpi_allreduce(norm,normmax,1,MPI_REAL_PRECISION,MPI_MAX,MPI_COMM_WORLD,mpierr)
#else /* WITH_MPI */
normmax = norm
#endif /* WITH_MPI */
#endif /* REALCASE == 1 */
if (myid .eq. 0) then
print *," Maximum error of result: ", normmax
endif
......
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