Commit 68d112f5 authored by Pavel Kus's avatar Pavel Kus Committed by Andreas Marek

adding scalapack test for part of the eigenvectors

parent 7dd73a2c
......@@ -14,6 +14,7 @@ solver_flag = {
"1stage" : "-DTEST_SOLVER_1STAGE",
"2stage" : "-DTEST_SOLVER_2STAGE",
"scalapack_all" : "-DTEST_SCALAPACK_ALL",
"scalapack_part" : "-DTEST_SCALAPACK_PART",
}
gpu_flag = {
0 : "-DTEST_GPU=0",
......@@ -50,7 +51,7 @@ for m, g, t, p, d, s, l in product(
if(m == "analytic" and (g == 1 or t != "eigenvectors")):
continue
if(s == "scalapack_all" and (g == 1 or t != "eigenvectors" or m != "analytic")):
if(s in ["scalapack_all", "scalapack_part"] and (g == 1 or t != "eigenvectors" or m != "analytic")):
continue
if (t == "solve_tridiagonal" and (s == "2stage" or d == "complex")):
......@@ -80,7 +81,7 @@ for m, g, t, p, d, s, l in product(
print("if WITH_MPI")
endifs += 1
if (s == "scalapack_all"):
if (s in ["scalapack_all", "scalapack_part"]):
print("if WITH_SCALAPACK_TESTS")
endifs += 1
......
......@@ -56,8 +56,8 @@ error: define exactly one of TEST_REAL or TEST_COMPLEX
error: define exactly one of TEST_SINGLE or TEST_DOUBLE
#endif
#if !(defined(TEST_SOLVER_1STAGE) ^ defined(TEST_SOLVER_2STAGE) ^ defined(TEST_SCALAPACK_ALL))
error: define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE or TEST_SCALAPACK_ALL
#if !(defined(TEST_SOLVER_1STAGE) ^ defined(TEST_SOLVER_2STAGE) ^ defined(TEST_SCALAPACK_ALL) ^ defined(TEST_SCALAPACK_PART))
error: define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE or TEST_SCALAPACK_ALL or TEST_SCALAPACK_PART
#endif
#ifdef TEST_SOLVER_1STAGE
......@@ -434,6 +434,8 @@ program test
call e%timer_start("e%eigenvectors()")
#ifdef TEST_SCALAPACK_ALL
call solve_scalapack_all(na, a, sc_desc, ev, z)
#elif TEST_SCALAPACK_PART
call solve_scalapack_part(na, a, sc_desc, nev, ev, z)
#else
call e%eigenvectors(a, ev, z, error)
#endif
......
......@@ -41,6 +41,7 @@
! any derivatives of ELPA under the same license that we chose for
! the original distribution, the GNU Lesser General Public License.
#include "../Fortran/assert.h"
#include "config-f90.h"
module test_scalapack
......@@ -57,6 +58,17 @@ module test_scalapack
#endif
end interface
interface solve_scalapack_part
module procedure solve_pdsyevr
module procedure solve_pzheevr
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure solve_pssyevr
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure solve_pcheevr
#endif
end interface
contains
#define COMPLEXCASE 1
......
......@@ -41,7 +41,7 @@
! any derivatives of ELPA under the same license that we chose for
! the original distribution, the GNU Lesser General Public License.
! compute all eigenvectores
subroutine solve_p&
&BLAS_CHAR_AND_SY_OR_HE&
&evd(na, a, sc_desc, ev, z)
......@@ -96,3 +96,66 @@
deallocate(iwork, work, rwork)
end subroutine
! compute part of eigenvectores
subroutine solve_p&
&BLAS_CHAR_AND_SY_OR_HE&
&evr(na, a, sc_desc, nev, ev, z)
implicit none
#include "../../src/general/precision_kinds.F90"
integer(kind=ik), intent(in) :: na, nev
MATH_DATATYPE(kind=rck), intent(in) :: a(:,:)
MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:)
real(kind=rk), intent(inout) :: ev(:)
integer(kind=ik), intent(in) :: sc_desc(:)
integer(kind=ik) :: info, lwork, liwork, lrwork
MATH_DATATYPE(kind=rck), allocatable :: work(:)
real(kind=rk), allocatable :: rwork(:)
integer, allocatable :: iwork(:)
integer(kind=ik) :: comp_eigenval, comp_eigenvec, smallest_ev_idx, largest_ev_idx
allocate(work(1), iwork(1), rwork(1))
smallest_ev_idx = 1
largest_ev_idx = nev
! query for required workspace
#ifdef REALCASE
call p&
&BLAS_CHAR&
&syevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, -1, iwork, -1, info)
#endif
#ifdef COMPLEXCASE
call p&
&BLAS_CHAR&
&heevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, -1, rwork, -1, iwork, -1, info)
#endif
! write(*,*) "computed sizes", lwork, liwork, "required sizes ", work(1), iwork(1)
lwork = work(1)
liwork = iwork(1)
deallocate(work, iwork)
allocate(work(lwork), stat = info)
allocate(iwork(liwork), stat = info)
#ifdef COMPLEXCASE
lrwork = rwork(1)
deallocate(rwork)
allocate(rwork(lrwork), stat = info)
#endif
! the actuall call to the method
#ifdef REALCASE
call p&
&BLAS_CHAR&
&syevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, lwork, iwork, liwork, info)
#endif
#ifdef COMPLEXCASE
call p&
&BLAS_CHAR&
&heevr('V', 'I', 'U', na, a, 1, 1, sc_desc, 0.0_rk, 0.0_rk, smallest_ev_idx, largest_ev_idx, &
comp_eigenval, comp_eigenvec, ev, z, 1, 1, sc_desc, work, lwork, rwork, lrwork, iwork, liwork, info)
#endif
assert(comp_eigenval == nev)
assert(comp_eigenvec == nev)
deallocate(iwork, work, rwork)
end subroutine
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