Commit 098bd541 authored by Pavel Kus's avatar Pavel Kus

progress in generalized evp

parent 2b5d44a3
......@@ -83,6 +83,7 @@ endif
EXTRA_libelpa@SUFFIX@_private_la_DEPENDENCIES = \
src/elpa1/elpa_reduce_add_vectors.F90 \
src/elpa1/elpa_transpose_vectors.F90 \
src/elpa_impl_template.F90 \
src/elpa1/elpa1_compute_template.F90 \
src/elpa2/elpa2_compute_real_template.F90 \
src/elpa2/elpa2_compute_complex_template.F90 \
......@@ -582,6 +583,7 @@ EXTRA_DIST = \
manual_cpp \
nvcc_wrap \
src/GPU/cuUtils_template.cu \
src/elpa_impl_template.F90 \
src/elpa1/elpa1_compute_template.F90 \
src/elpa1/elpa1_merge_systems_real_template.F90 \
src/elpa1/elpa1_solve_tridi_real_template.F90 \
......
......@@ -100,6 +100,7 @@ test_type_flag = {
"cholesky" : "-DTEST_CHOLESKY",
"hermitian_multiply" : "-DTEST_HERMITIAN_MULTIPLY",
"qr" : "-DTEST_QR_DECOMPOSITION",
"generalized" : "-DTEST_GENERALIZED_EIGENPROBLEM",
}
layout_flag = {
......
......@@ -114,6 +114,12 @@ module elpa_api
elpa_eigenvalues_dc, &
elpa_eigenvalues_fc
generic, public :: generalized_eigenvectors => & !< method eigenvectors for solving the full generalized eigenvalue problem
elpa_generalized_eigenvectors_d, & !< the eigenvalues and (parts of) the eigenvectors are computed
elpa_generalized_eigenvectors_f, & !< for symmetric real valued / hermitian complex valued matrices
elpa_generalized_eigenvectors_dc, &
elpa_generalized_eigenvectors_fc
generic, public :: hermitian_multiply => & !< method for a "hermitian" multiplication of matrices a and b
elpa_hermitian_multiply_d, & !< for real valued matrices: a**T * b
elpa_hermitian_multiply_dc, & !< for complex valued matrices a**H * b
......@@ -154,6 +160,11 @@ module elpa_api
procedure(elpa_eigenvalues_dc_i), deferred, public :: elpa_eigenvalues_dc
procedure(elpa_eigenvalues_fc_i), deferred, public :: elpa_eigenvalues_fc
procedure(elpa_generalized_eigenvectors_d_i), deferred, public :: elpa_generalized_eigenvectors_d
procedure(elpa_generalized_eigenvectors_f_i), deferred, public :: elpa_generalized_eigenvectors_f
procedure(elpa_generalized_eigenvectors_dc_i), deferred, public :: elpa_generalized_eigenvectors_dc
procedure(elpa_generalized_eigenvectors_fc_i), deferred, public :: elpa_generalized_eigenvectors_fc
procedure(elpa_hermitian_multiply_d_i), deferred, public :: elpa_hermitian_multiply_d
procedure(elpa_hermitian_multiply_f_i), deferred, public :: elpa_hermitian_multiply_f
procedure(elpa_hermitian_multiply_dc_i), deferred, public :: elpa_hermitian_multiply_dc
......@@ -661,6 +672,156 @@ module elpa_api
end subroutine
end interface
!> \brief abstract definition of interface to solve double real generalized eigenvalue problem
!>
!> The dimensions of the matrix a and b (locally ditributed and global), the block-cyclic distribution
!> blocksize, the number of eigenvectors
!> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
!> with the class method "setup"
!>
!> It is possible to change the behaviour of the method by setting tunable parameters with the
!> class method "set"
!> Parameters
!> \details
!> \param self class(elpa_t), the ELPA object
!> \param a double real matrix a: defines the problem to solve
!> \param b double real matrix b: defines the problem to solve
!> \param ev double real: on output stores the eigenvalues
!> \param q double real matrix q: on output stores the eigenvalues
!> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface
subroutine elpa_generalized_eigenvectors_d_i(self, a, b, ev, q, sc_desc, error)
use iso_c_binding
use elpa_constants
import elpa_t
implicit none
class(elpa_t) :: self
#ifdef USE_ASSUMED_SIZE
real(kind=c_double) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *)
#else
real(kind=c_double) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), &
q(self%local_nrows, self%local_ncols)
#endif
real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
end subroutine
end interface
!> \brief abstract definition of interface to solve single real generalized eigenvalue problem
!>
!> The dimensions of the matrix a and b(locally ditributed and global), the block-cyclic distribution
!> blocksize, the number of eigenvectors
!> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
!> with the class method "setup"
!>
!> It is possible to change the behaviour of the method by setting tunable parameters with the
!> class method "set"
!> Parameters
!> \details
!> \param self class(elpa_t), the ELPA object
!> \param a single real matrix a: defines the problem to solve
!> \param b single real matrix b: defines the problem to solve
!> \param ev single real: on output stores the eigenvalues
!> \param q single real matrix q: on output stores the eigenvalues
!> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface
subroutine elpa_generalized_eigenvectors_f_i(self, a, b, ev, q, sc_desc, error)
use iso_c_binding
use elpa_constants
import elpa_t
implicit none
class(elpa_t) :: self
#ifdef USE_ASSUMED_SIZE
real(kind=c_float) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *)
#else
real(kind=c_float) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), &
q(self%local_nrows, self%local_ncols)
#endif
real(kind=c_float) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
end subroutine
end interface
!> \brief abstract definition of interface to solve double complex generalized eigenvalue problem
!>
!> The dimensions of the matrix a and b(locally ditributed and global), the block-cyclic distribution
!> blocksize, the number of eigenvectors
!> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
!> with the class method "setup"
!>
!> It is possible to change the behaviour of the method by setting tunable parameters with the
!> class method "set"
!> Parameters
!> \details
!> \param self class(elpa_t), the ELPA object
!> \param a double complex matrix a: defines the problem to solve
!> \param b double complex matrix b: defines the problem to solve
!> \param ev double real: on output stores the eigenvalues
!> \param q double complex matrix q: on output stores the eigenvalues
!> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface
subroutine elpa_generalized_eigenvectors_dc_i(self, a, b, ev, q, sc_desc, error)
use iso_c_binding
use elpa_constants
import elpa_t
implicit none
class(elpa_t) :: self
#ifdef USE_ASSUMED_SIZE
complex(kind=c_double_complex) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *)
#else
complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), &
q(self%local_nrows, self%local_ncols)
#endif
real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
end subroutine
end interface
!> \brief abstract definition of interface to solve single complex generalized eigenvalue problem
!>
!> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution
!> blocksize, the number of eigenvectors
!> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE
!> with the class method "setup"
!>
!> It is possible to change the behaviour of the method by setting tunable parameters with the
!> class method "set"
!> Parameters
!> \details
!> \param self class(elpa_t), the ELPA object
!> \param a single complex matrix a: defines the problem to solve
!> \param b single complex matrix b: defines the problem to solve
!> \param ev single real: on output stores the eigenvalues
!> \param q single complex matrix q: on output stores the eigenvalues
!> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface
subroutine elpa_generalized_eigenvectors_fc_i(self, a, b, ev, q, sc_desc, error)
use iso_c_binding
use elpa_constants
import elpa_t
implicit none
class(elpa_t) :: self
#ifdef USE_ASSUMED_SIZE
complex(kind=c_float_complex) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *)
#else
complex(kind=c_float_complex) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), &
q(self%local_nrows, self%local_ncols)
#endif
real(kind=c_float) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
end subroutine
end interface
!> \brief abstract definition of interface to compute C : = A**T * B for double real matrices
!> where A is a square matrix (self%a,self%na) which is optionally upper or lower triangular
!> B is a (self%na,ncb) matrix
......
......@@ -49,5 +49,8 @@ module elpa_constants
use, intrinsic :: iso_c_binding, only : C_INT
implicit none
public
integer(kind=C_INT), parameter :: SC_DESC_LEN = 9
#include "src/fortran_constants.F90"
end module
This diff is collapsed.
subroutine elpa_transform_generalized_&
&ELPA_IMPL_SUFFIX&
&(self, a, b, sc_desc, error)
implicit none
#include "general/precision_kinds.F90"
class(elpa_impl_t) :: self
#ifdef USE_ASSUMED_SIZE
MATH_DATATYPE(kind=rck) :: a(self%local_nrows, *), b(self%local_nrows, *)
#else
MATH_DATATYPE(kind=rck) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols)
#endif
integer :: error
integer :: sc_desc(9)
! 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)
! TODO: why I cannot use self%elpa ??
! B = U^T*U, B<-U
call self%elpa_cholesky_&
&ELPA_IMPL_SUFFIX&
&(b, error)
if(error .NE. ELPA_OK) return
! B <- inv(U)
call self%elpa_invert_trm_&
&ELPA_IMPL_SUFFIX&
&(b, error)
if(error .NE. ELPA_OK) return
! ! tmp <- inv(U^T) * A
! call self%elpa_hermitian_multiply_&
! &ELPA_IMPL_SUFFIX&
! &('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 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)
! A <= inv(U)^T * A * inv(U)
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
end subroutine
......@@ -43,6 +43,7 @@
#undef scal_PRECISION_NRM2
#undef scal_PRECISION_LASET
#undef PRECISION_SUFFIX
#undef ELPA_IMPL_SUFFIX
#undef MPI_REAL_PRECISION
#undef MPI_MATH_DATATYPE_PRECISION
......@@ -60,6 +61,7 @@
#define PRECISION double
#define PRECISION_STR 'double'
#define PRECISION_SUFFIX "_double"
#define ELPA_IMPL_SUFFIX d
#define REAL_DATATYPE rk8
#define BLAS_CHAR D
#define BLAS_CHAR_AND_SY_OR_HE DSY
......@@ -112,6 +114,7 @@
#define PRECISION single
#define PRECISION_STR 'single'
#define PRECISION_SUFFIX "_single"
#define ELPA_IMPL_SUFFIX f
#define REAL_DATATYPE rk4
#define BLAS_CHAR S
#define BLAS_CHAR_AND_SY_OR_HE SSY
......@@ -206,6 +209,7 @@
#undef scal_PRECISION_DOTC
#undef scal_PRECISION_LASET
#undef PRECISION_SUFFIX
#undef ELPA_IMPL_SUFFIX
#undef MPI_COMPLEX_PRECISION
#undef MPI_MATH_DATATYPE_PRECISION
#undef MPI_MATH_DATATYPE_PRECISION_EXPL
......@@ -229,6 +233,7 @@
#define PRECISION double
#define PRECISION_STR 'double'
#define PRECISION_SUFFIX "_double"
#define ELPA_IMPL_SUFFIX dc
#define COMPLEX_DATATYPE CK8
#define BLAS_CHAR Z
#define BLAS_CHAR_AND_SY_OR_HE ZHE
......@@ -285,6 +290,7 @@
#define PRECISION single
#define PRECISION_STR 'single'
#define PRECISION_SUFFIX "_single"
#define ELPA_IMPL_SUFFIX fc
#define COMPLEX_DATATYPE CK4
#define BLAS_CHAR C
#define BLAS_CHAR_AND_SY_OR_HE CHE
......
......@@ -87,12 +87,22 @@ error: define either TEST_ALL_KERNELS or a valid TEST_KERNEL
# define REALCASE 1
# undef COMPLEXCASE
# define MATH_DATATYPE real
#define KERNEL_KEY "real_kernel"
# define KERNEL_KEY "real_kernel"
# ifdef TEST_SINGLE
# define BLAS_CHAR S
# else
# define BLAS_CHAR D
# endif
#else
# define COMPLEXCASE 1
# undef REALCASE
# define MATH_DATATYPE complex
#define KERNEL_KEY "complex_kernel"
# define KERNEL_KEY "complex_kernel"
# ifdef TEST_SINGLE
# define BLAS_CHAR C
# else
# define BLAS_CHAR Z
# endif
#endif
#include "assert.h"
......@@ -138,6 +148,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(:,:)
......@@ -146,14 +159,8 @@ program test
logical :: check_all_evals
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION) || defined(TEST_HERMITIAN_MULTIPLY)
EV_TYPE, allocatable :: d(:), sd(:), ds(:), sds(:)
EV_TYPE :: diagonalELement, subdiagonalElement
#endif
#if defined(TEST_CHOLESKY)
MATRIX_TYPE, allocatable :: d(:), sd(:), ds(:), sds(:)
MATRIX_TYPE :: diagonalELement, subdiagonalElement
#endif
integer :: error, status
......@@ -274,6 +281,7 @@ program test
#endif
call e%set("timings",1)
call e%set("debug", 1)
assert_elpa_ok(e%setup())
......@@ -305,6 +313,14 @@ 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))
! todo: only temporarily, before we start using random SPD matrix for B
allocate(d (na), ds(na))
allocate(sd (na), sds(na))
#endif
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION) || defined(TEST_CHOLESKY)
allocate(d (na), ds(na))
allocate(sd (na), sds(na))
......@@ -312,18 +328,28 @@ program test
#endif
! prepare matrices
call e%timer_start("e%prepare_matrices()")
call e%timer_start("prepare_matrices")
a(:,:) = 0.0
z(:,:) = 0.0
ev(:) = 0.0
#if defined(TEST_EIGENVECTORS) || defined(TEST_HERMITIAN_MULTIPLY) || defined(TEST_QR_DECOMPOSITION)
#if defined(TEST_EIGENVECTORS) || defined(TEST_HERMITIAN_MULTIPLY) || defined(TEST_QR_DECOMPOSITION) || defined(TEST_GENERALIZED_EIGENPROBLEM)
#ifdef TEST_MATRIX_ANALYTIC
call prepare_matrix_analytic(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
as(:,:) = a
#else
#else /* TEST_MATRIX_ANALYTIC */
if (nev .ge. 1) then
call prepare_matrix_random(na, myid, sc_desc, a, z, as)
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
!call prepare_matrix_random(na, myid, sc_desc, b, z, bs)
! TODO create random SPD matrix
!diagonalElement = (2.546_rk, 0.0_rk)
diagonalElement = ONE
subdiagonalElement = (0.0_rk, 0.0_rk)
call prepare_matrix_toeplitz(na, diagonalElement, subdiagonalElement, &
d, sd, ds, sds, b, bs, nblk, np_rows, &
np_cols, my_prow, my_pcol)
#endif
else
! zero eigenvectors and not analytic test => toeplitz matrix
diagonalElement = 0.45_rk
......@@ -338,7 +364,7 @@ program test
c(:,:) = ONE
#endif /* TEST_HERMITIAN_MULTIPLY */
#endif /* TEST_MATRIX_ANALYTIC */
#endif /* (not) TEST_MATRIX_ANALYTIC */
#endif /* defined(TEST_EIGENVECTORS) || defined(TEST_HERMITIAN_MULTIPLY) || defined(TEST_QR_DECOMPOSITION) */
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL)
......@@ -356,10 +382,10 @@ program test
d, sd, ds, sds, a, as, nblk, np_rows, &
np_cols, my_prow, my_pcol)
#endif /* TEST_CHOLESKY */
call e%timer_stop("e%prepare_matrices()")
call e%timer_stop("prepare_matrices")
if (myid == 0) then
print *, ""
call e%print_times("e%prepare_matrices()")
call e%print_times("prepare_matrices")
print *, ""
endif
......@@ -392,20 +418,42 @@ program test
call e%timer_start(elpa_int_value_to_string(KERNEL_KEY, kernel))
#endif
!#if defined(TEST_GENERALIZED_EIGENPROBLEM)
! call e%timer_start("generalized_eigenproblem")
! call e%timer_start("e%cholesky()")
! call e%cholesky(b, error)
! assert_elpa_ok(error)
! call e%timer_stop("e%cholesky()")
! call e%timer_start("e%invert_triangular")
! call e%invert_triangular(b, error)
! assert_elpa_ok(error)
! call e%timer_stop("e%invert_triangular")
!#ifdef WITH_MPI
!
!#endif
!#endif
! The actual solve step
#if defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION)
#if defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION) || defined(TEST_GENERALIZED_EIGENPROBLEM)
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)
check_all_evals = .false. ! scalapack does not compute all eigenvectors
#elif TEST_GENERALIZED_EIGENPROBLEM
call e%generalized_eigenvectors(a, b, ev, z, sc_desc, error)
#else
call e%eigenvectors(a, ev, z, error)
#endif
call e%timer_stop("e%eigenvectors()")
#endif /* TEST_EIGENVECTORS || defined(TEST_QR_DECOMPOSITION) */
!#if defined(TEST_GENERALIZED_EIGENPROBLEM)
! ! todo...
! call e%timer_stop("generalized_eigenproblem")
!#endif
#ifdef TEST_EIGENVALUES
call e%timer_start("e%eigenvalues()")
call e%eigenvalues(a, ev, error)
......@@ -442,7 +490,7 @@ program test
call e%print_times(elpa_int_value_to_string(KERNEL_KEY, kernel))
#else /* TEST_ALL_KERNELS */
#if defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION)
#if defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION) || defined(TEST_GENERALIZED_EIGENPROBLEM)
call e%print_times("e%eigenvectors()")
#endif
#ifdef TEST_EIGENVALUES
......@@ -461,8 +509,8 @@ program test
endif
! check the results
call e%timer_start("e%check_correctness()")
#if defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION)
call e%timer_start("check_correctness")
#if defined(TEST_EIGENVECTORS) || defined(TEST_QR_DECOMPOSITION) || defined(TEST_GENERALIZED_EIGENPROBLEM)
#ifdef TEST_MATRIX_ANALYTIC
status = check_correctness_analytic(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals)
#else
......@@ -504,7 +552,7 @@ program test
status = check_correctness_hermitian_multiply(na, a, b, c, na_rows, sc_desc, myid )
call check_status(status, myid)
#endif
call e%timer_stop("e%check_correctness()")
call e%timer_stop("check_correctness")
if (myid == 0) then
print *, ""
......@@ -521,7 +569,7 @@ program test
if (myid == 0) then
print *, ""
call e%print_times("e%check_correctness()")
call e%print_times("check_correctness")
print *, ""
endif
......@@ -543,6 +591,14 @@ program test
deallocate(ev_analytic)
#endif
#ifdef TEST_GENERALIZED_EIGENPROBLEM
deallocate(b)
deallocate(bs)
! todo: only temporarily, before we start using random SPD matrix for B
deallocate(d, ds)
deallocate(sd, sds)
#endif
#ifdef TEST_ALL_LAYOUTS
end do ! factors
end do ! layouts
......
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