Commit 316c3a39 authored by Pavel Kus's avatar Pavel Kus

re-enabled generalized eigenproblem test

* tested with elpa1 only (should be sufficient)
* not using Hermitian multiply at the moment
* requires scalapack
parent 1f2e18c5
......@@ -44,6 +44,7 @@ test_type_flag = {
"solve_tridiagonal": "-DTEST_SOLVE_TRIDIAGONAL",
"cholesky": "-DTEST_CHOLESKY",
"hermitian_multiply": "-DTEST_HERMITIAN_MULTIPLY",
"generalized" : "-DTEST_GENERALIZED_EIGENPROBLEM",
}
layout_flag = {
......@@ -89,6 +90,10 @@ for lang, m, g, q, t, p, d, s, lay in product(sorted(language_flag.keys()),
if (t == "solve_tridiagonal" and (s != "1stage" or d != "real" or m != "toeplitz")):
continue
# solve generalized only for random matrix in 1stage
if (t == "generalized" and (m != "random" or s == "2stage")):
continue
# cholesky tests only 1stage and teoplitz or random matrix
if (t == "cholesky" and ((not (m == "toeplitz" or m == "random")) or s == "2stage")):
continue
......
......@@ -16,6 +16,7 @@
! local helper array. TODO: do we want it this way? (do we need it? )
MATH_DATATYPE(kind=rck) :: tmp(self%local_nrows, self%local_ncols)
call self%timer_start("transform_generalized()")
! TODO: why I cannot use self%elpa ??
! B = U^T*U, B<-U
call self%elpa_cholesky_&
......@@ -33,26 +34,36 @@
! &('U','F', self%na, b, a, self%local_nrows, self%local_ncols, tmp, &
! self%local_nrows, self%local_ncols, error)
! if(error .NE. ELPA_OK) return
#ifdef WITH_MPI
! A <= inv(U)^T * A
call self%timer_start("scalapack multiply inv(U)^T * A")
#ifdef WITH_MPI
call p&
&BLAS_CHAR&
&trmm("L", "U", BLAS_TRANS_OR_CONJ, "N", self%na, self%na, &
ONE, b, 1, 1, sc_desc, a, 1, 1, sc_desc)
#else
call BLAS_CHAR&
&trmm("L", "U", BLAS_TRANS_OR_CONJ, "N", self%na, self%na, &
ONE, b, self%na, a, self%na)
#endif
call self%timer_stop("scalapack multiply inv(U)^T * A")
! A <= inv(U)^T * A * inv(U)
call self%timer_start("scalapack multiply A * inv(U)")
#ifdef WITH_MPI
call p&
&BLAS_CHAR&
&trmm("R", "U", "N", "N", self%na, self%na, &
ONE, b, 1, 1, sc_desc, a, 1, 1, sc_desc)
#else
call BLAS_CHAR&
&trmm("L", "U", BLAS_TRANS_OR_CONJ, "N", self%na, self%na, &
ONE, b, self%na, a, self%na)
call BLAS_CHAR&
&trmm("R", "U", "N", "N", self%na, self%na, &
ONE, b, self%na, a, self%na)
#endif
call self%timer_stop("scalapack multiply A * inv(U)")
call self%timer_stop("transform_generalized()")
end subroutine
......@@ -73,7 +84,10 @@
! local helper array. TODO: do we want it this way? (do we need it? )
MATH_DATATYPE(kind=rck) :: tmp(self%local_nrows, self%local_ncols)
call self%timer_start("transform_back_generalized()")
!todo: part of eigenvectors only
call self%timer_start("scalapack multiply inv(U) * Q")
#ifdef WITH_MPI
! Q <= inv(U) * Q
call p&
......@@ -85,6 +99,9 @@
&trmm("L", "U", "N", "N", self%na, self%na, &
ONE, b, self%na, q, self%na)
#endif
call self%timer_stop("scalapack multiply inv(U) * Q")
call self%timer_stop("transform_back_generalized()")
end subroutine
#endif
......
......@@ -137,6 +137,9 @@ program test
MATRIX_TYPE, allocatable :: a(:,:), as(:,:)
#if defined(TEST_HERMITIAN_MULTIPLY)
MATRIX_TYPE, allocatable :: b(:,:), c(:,:)
#endif
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
MATRIX_TYPE, allocatable :: b(:,:), bs(:,:)
#endif
! eigenvectors
MATRIX_TYPE, allocatable :: z(:,:)
......@@ -163,7 +166,8 @@ program test
#endif
integer :: kernel
character(len=1) :: layout
logical :: do_test_numeric_residual, do_test_analytic_eigenvalues, &
logical :: do_test_numeric_residual, do_test_numeric_residual_generalized, &
do_test_analytic_eigenvalues, &
do_test_analytic_eigenvalues_eigenvectors, &
do_test_frank_eigenvalues, &
do_test_toeplitz_eigenvalues, do_test_cholesky, &
......@@ -182,6 +186,7 @@ program test
do_test_numeric_residual = .false.
do_test_numeric_residual_generalized = .false.
do_test_analytic_eigenvalues = .false.
do_test_analytic_eigenvalues_eigenvectors = .false.
do_test_frank_eigenvalues = .false.
......@@ -276,6 +281,11 @@ program test
allocate(c (na_rows,na_cols))
#endif
#ifdef TEST_GENERALIZED_EIGENPROBLEM
allocate(b (na_rows,na_cols))
allocate(bs (na_rows,na_cols))
#endif
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
allocate(d (na), ds(na))
allocate(sd (na), sds(na))
......@@ -327,6 +337,18 @@ program test
do_test_toeplitz_eigenvalues = .false.
#endif /* TEST_MATRIX_RANDOM and TEST_CHOLESKY */
#if defined(TEST_MATRIX_RANDOM) && defined(TEST_GENERALIZED_EIGENPROBLEM)
! call prepare_matrix_random(na, myid, sc_desc, a, z, as)
call prepare_matrix_random_spd(na, myid, sc_desc, b, z, bs, &
nblk, np_rows, np_cols, my_prow, my_pcol)
do_test_analytic_eigenvalues = .false.
do_test_analytic_eigenvalues_eigenvectors = .false.
do_test_frank_eigenvalues = .false.
do_test_toeplitz_eigenvalues = .false.
do_test_numeric_residual = .false.
do_test_numeric_residual_generalized = .true.
#endif /* TEST_MATRIX_RANDOM and TEST_GENERALIZED_EIGENPROBLEM */
#if defined(TEST_MATRIX_RANDOM) && (defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVALUES))
#error "Random matrix is not allowed in this configuration"
#endif
......@@ -601,6 +623,12 @@ program test
call e%timer_stop("e%hermitian_multiply()")
#endif
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
call e%timer_start("e%generalized_eigenvectors()")
call e%generalized_eigenvectors(a, b, ev, z, sc_desc, error)
call e%timer_stop("e%generalized_eigenvectors()")
#endif
assert_elpa_ok(error)
#ifdef TEST_ALL_KERNELS
......@@ -631,6 +659,9 @@ program test
#ifdef TEST_HERMITIAN_MULTIPLY
call e%print_times("e%hermitian_multiply()")
#endif
#ifdef TEST_GENERALIZED_EIGENPROBLEM
call e%print_times("e%generalized_eigenvectors()")
#endif
#endif /* TEST_ALL_KERNELS */
endif
......@@ -643,6 +674,7 @@ program test
status = check_correctness_analytic(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals, .true.)
call check_status(status, myid)
endif
if(do_test_numeric_residual) then
status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
call check_status(status, myid)
......@@ -673,6 +705,14 @@ program test
endif
#endif
#ifdef TEST_GENERALIZED_EIGENPROBLEM
if(do_test_numeric_residual_generalized) then
status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, &
my_pcol, bs)
call check_status(status, myid)
endif
#endif
if (myid == 0) then
print *, ""
endif
......@@ -700,6 +740,9 @@ program test
deallocate(d, ds)
deallocate(sd, sds)
#endif
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
deallocate(b, bs)
#endif
#ifdef TEST_ALL_LAYOUTS
end do ! factors
......
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