Commit e93f1c06 authored by Pavel Kus's avatar Pavel Kus

generalized eigenproblem reuse of B

added possibility to re-use the already decomposed matrix B in the
generalized problem
Added respective test
parent 50d499f0
...@@ -45,6 +45,7 @@ test_type_flag = { ...@@ -45,6 +45,7 @@ test_type_flag = {
"cholesky": "-DTEST_CHOLESKY", "cholesky": "-DTEST_CHOLESKY",
"hermitian_multiply": "-DTEST_HERMITIAN_MULTIPLY", "hermitian_multiply": "-DTEST_HERMITIAN_MULTIPLY",
"generalized" : "-DTEST_GENERALIZED_EIGENPROBLEM", "generalized" : "-DTEST_GENERALIZED_EIGENPROBLEM",
"generalized_decomp": "-DTEST_GENERALIZED_DECOMP_EIGENPROBLEM",
} }
layout_flag = { layout_flag = {
...@@ -94,6 +95,11 @@ for lang, m, g, q, t, p, d, s, lay in product(sorted(language_flag.keys()), ...@@ -94,6 +95,11 @@ for lang, m, g, q, t, p, d, s, lay in product(sorted(language_flag.keys()),
if (t == "generalized" and (m != "random" or s == "2stage")): if (t == "generalized" and (m != "random" or s == "2stage")):
continue continue
# solve generalized already decomposed only for random matrix in 1stage
# maybe this test should be further restricted, maybe not so important...
if (t == "generalized_decomp" and (m != "random" or s == "2stage")):
continue
# cholesky tests only 1stage and teoplitz or random matrix # cholesky tests only 1stage and teoplitz or random matrix
if (t == "cholesky" and ((not (m == "toeplitz" or m == "random")) or s == "2stage")): if (t == "cholesky" and ((not (m == "toeplitz" or m == "random")) or s == "2stage")):
continue continue
......
...@@ -821,9 +821,10 @@ module elpa_api ...@@ -821,9 +821,10 @@ module elpa_api
!> \param b double real matrix b: 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 ev double real: on output stores the eigenvalues
!> \param q double real matrix q: on output stores the eigenvalues !> \param q double real matrix q: on output stores the eigenvalues
!> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)?
!> \result error integer, optional : error code, which can be queried with elpa_strerr !> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface abstract interface
subroutine elpa_generalized_eigenvectors_d_i(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_d_i(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use iso_c_binding use iso_c_binding
use elpa_constants use elpa_constants
import elpa_t import elpa_t
...@@ -837,7 +838,7 @@ module elpa_api ...@@ -837,7 +838,7 @@ module elpa_api
#endif #endif
real(kind=c_double) :: ev(self%na) real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
logical :: is_already_decomposed
integer, optional :: error integer, optional :: error
end subroutine end subroutine
end interface end interface
...@@ -858,9 +859,10 @@ module elpa_api ...@@ -858,9 +859,10 @@ module elpa_api
!> \param b single real matrix b: 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 ev single real: on output stores the eigenvalues
!> \param q single real matrix q: on output stores the eigenvalues !> \param q single real matrix q: on output stores the eigenvalues
!> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)?
!> \result error integer, optional : error code, which can be queried with elpa_strerr !> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface abstract interface
subroutine elpa_generalized_eigenvectors_f_i(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_f_i(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use iso_c_binding use iso_c_binding
use elpa_constants use elpa_constants
import elpa_t import elpa_t
...@@ -874,6 +876,7 @@ module elpa_api ...@@ -874,6 +876,7 @@ module elpa_api
#endif #endif
real(kind=c_float) :: ev(self%na) real(kind=c_float) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
logical :: is_already_decomposed
integer, optional :: error integer, optional :: error
end subroutine end subroutine
...@@ -895,9 +898,10 @@ module elpa_api ...@@ -895,9 +898,10 @@ module elpa_api
!> \param b double complex matrix b: 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 ev double real: on output stores the eigenvalues
!> \param q double complex matrix q: on output stores the eigenvalues !> \param q double complex matrix q: on output stores the eigenvalues
!> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)?
!> \result error integer, optional : error code, which can be queried with elpa_strerr !> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface abstract interface
subroutine elpa_generalized_eigenvectors_dc_i(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_dc_i(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use iso_c_binding use iso_c_binding
use elpa_constants use elpa_constants
import elpa_t import elpa_t
...@@ -912,6 +916,7 @@ module elpa_api ...@@ -912,6 +916,7 @@ module elpa_api
#endif #endif
real(kind=c_double) :: ev(self%na) real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
logical :: is_already_decomposed
integer, optional :: error integer, optional :: error
end subroutine end subroutine
...@@ -933,9 +938,10 @@ module elpa_api ...@@ -933,9 +938,10 @@ module elpa_api
!> \param b single complex matrix b: 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 ev single real: on output stores the eigenvalues
!> \param q single complex matrix q: on output stores the eigenvalues !> \param q single complex matrix q: on output stores the eigenvalues
!> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)?
!> \result error integer, optional : error code, which can be queried with elpa_strerr !> \result error integer, optional : error code, which can be queried with elpa_strerr
abstract interface abstract interface
subroutine elpa_generalized_eigenvectors_fc_i(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_fc_i(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use iso_c_binding use iso_c_binding
use elpa_constants use elpa_constants
import elpa_t import elpa_t
...@@ -949,6 +955,7 @@ module elpa_api ...@@ -949,6 +955,7 @@ module elpa_api
#endif #endif
real(kind=c_float) :: ev(self%na) real(kind=c_float) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
logical :: is_already_decomposed
integer, optional :: error integer, optional :: error
end subroutine end subroutine
......
...@@ -1511,7 +1511,7 @@ module elpa_impl ...@@ -1511,7 +1511,7 @@ module elpa_impl
!> even if only a part of the eigenvalues is needed. !> even if only a part of the eigenvalues is needed.
!> !>
!> \param error integer, optional: returns an error code, which can be queried with elpa_strerr !> \param error integer, optional: returns an error code, which can be queried with elpa_strerr
subroutine elpa_generalized_eigenvectors_d(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_d(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use elpa2_impl use elpa2_impl
use elpa1_impl use elpa1_impl
use elpa_utilities, only : error_unit use elpa_utilities, only : error_unit
...@@ -1526,13 +1526,14 @@ module elpa_impl ...@@ -1526,13 +1526,14 @@ module elpa_impl
#endif #endif
real(kind=c_double) :: ev(self%na) real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
logical :: is_already_decomposed
integer, optional :: error integer, optional :: error
integer :: error_l integer :: error_l
integer(kind=c_int) :: solver integer(kind=c_int) :: solver
logical :: success_l logical :: success_l
call self%elpa_transform_generalized_d(a, b, sc_desc, error_l) call self%elpa_transform_generalized_d(a, b, sc_desc, is_already_decomposed, error_l)
if (present(error)) then if (present(error)) then
error = error_l error = error_l
else if (error_l .ne. ELPA_OK) then else if (error_l .ne. ELPA_OK) then
...@@ -1585,7 +1586,7 @@ module elpa_impl ...@@ -1585,7 +1586,7 @@ module elpa_impl
call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols]) call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN]) call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN])
call elpa_generalized_eigenvectors_d(self, a, b, ev, q, sc_desc, error) call elpa_generalized_eigenvectors_d(self, a, b, ev, q, sc_desc, .false., error)
end subroutine end subroutine
...@@ -1614,7 +1615,7 @@ module elpa_impl ...@@ -1614,7 +1615,7 @@ module elpa_impl
!> even if only a part of the eigenvalues is needed. !> even if only a part of the eigenvalues is needed.
!> !>
!> \param error integer, optional: returns an error code, which can be queried with elpa_strerr !> \param error integer, optional: returns an error code, which can be queried with elpa_strerr
subroutine elpa_generalized_eigenvectors_f(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_f(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use elpa2_impl use elpa2_impl
use elpa1_impl use elpa1_impl
use elpa_utilities, only : error_unit use elpa_utilities, only : error_unit
...@@ -1630,12 +1631,13 @@ module elpa_impl ...@@ -1630,12 +1631,13 @@ module elpa_impl
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error integer, optional :: error
logical :: is_already_decomposed
integer :: error_l integer :: error_l
integer(kind=c_int) :: solver integer(kind=c_int) :: solver
#ifdef WANT_SINGLE_PRECISION_REAL #ifdef WANT_SINGLE_PRECISION_REAL
logical :: success_l logical :: success_l
call self%elpa_transform_generalized_f(a, b, sc_desc, error_l) call self%elpa_transform_generalized_f(a, b, sc_desc, is_already_decomposed, error_l)
if (present(error)) then if (present(error)) then
error = error_l error = error_l
else if (error_l .ne. ELPA_OK) then else if (error_l .ne. ELPA_OK) then
...@@ -1693,7 +1695,7 @@ module elpa_impl ...@@ -1693,7 +1695,7 @@ module elpa_impl
call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols]) call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN]) call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN])
call elpa_generalized_eigenvectors_f(self, a, b, ev, q, sc_desc, error) call elpa_generalized_eigenvectors_f(self, a, b, ev, q, sc_desc, .false., error)
end subroutine end subroutine
...@@ -1722,7 +1724,7 @@ module elpa_impl ...@@ -1722,7 +1724,7 @@ module elpa_impl
!> even if only a part of the eigenvalues is needed. !> even if only a part of the eigenvalues is needed.
!> !>
!> \param error integer, optional: returns an error code, which can be queried with elpa_strerr !> \param error integer, optional: returns an error code, which can be queried with elpa_strerr
subroutine elpa_generalized_eigenvectors_dc(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_dc(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use elpa2_impl use elpa2_impl
use elpa1_impl use elpa1_impl
use elpa_utilities, only : error_unit use elpa_utilities, only : error_unit
...@@ -1739,11 +1741,12 @@ module elpa_impl ...@@ -1739,11 +1741,12 @@ module elpa_impl
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error integer, optional :: error
logical :: is_already_decomposed
integer :: error_l integer :: error_l
integer(kind=c_int) :: solver integer(kind=c_int) :: solver
logical :: success_l logical :: success_l
call self%elpa_transform_generalized_dc(a, b, sc_desc, error_l) call self%elpa_transform_generalized_dc(a, b, sc_desc, is_already_decomposed, error_l)
if (present(error)) then if (present(error)) then
error = error_l error = error_l
else if (error_l .ne. ELPA_OK) then else if (error_l .ne. ELPA_OK) then
...@@ -1798,7 +1801,7 @@ module elpa_impl ...@@ -1798,7 +1801,7 @@ module elpa_impl
call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols]) call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN]) call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN])
call elpa_generalized_eigenvectors_dc(self, a, b, ev, q, sc_desc, error) call elpa_generalized_eigenvectors_dc(self, a, b, ev, q, sc_desc, .false., error)
end subroutine end subroutine
...@@ -1827,7 +1830,7 @@ module elpa_impl ...@@ -1827,7 +1830,7 @@ module elpa_impl
!> even if only a part of the eigenvalues is needed. !> even if only a part of the eigenvalues is needed.
!> !>
!> \param error integer, optional: returns an error code, which can be queried with elpa_strerr !> \param error integer, optional: returns an error code, which can be queried with elpa_strerr
subroutine elpa_generalized_eigenvectors_fc(self, a, b, ev, q, sc_desc, error) subroutine elpa_generalized_eigenvectors_fc(self, a, b, ev, q, sc_desc, is_already_decomposed, error)
use elpa2_impl use elpa2_impl
use elpa1_impl use elpa1_impl
use elpa_utilities, only : error_unit use elpa_utilities, only : error_unit
...@@ -1844,12 +1847,13 @@ module elpa_impl ...@@ -1844,12 +1847,13 @@ module elpa_impl
integer :: sc_desc(SC_DESC_LEN) integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error integer, optional :: error
logical :: is_already_decomposed
integer :: error_l integer :: error_l
integer(kind=c_int) :: solver integer(kind=c_int) :: solver
#ifdef WANT_SINGLE_PRECISION_COMPLEX #ifdef WANT_SINGLE_PRECISION_COMPLEX
logical :: success_l logical :: success_l
call self%elpa_transform_generalized_fc(a, b, sc_desc, error_l) call self%elpa_transform_generalized_fc(a, b, sc_desc, is_already_decomposed, error_l)
if (present(error)) then if (present(error)) then
error = error_l error = error_l
else if (error_l .ne. ELPA_OK) then else if (error_l .ne. ELPA_OK) then
...@@ -1908,7 +1912,7 @@ module elpa_impl ...@@ -1908,7 +1912,7 @@ module elpa_impl
call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols]) call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN]) call c_f_pointer(sc_desc_p, sc_desc, [SC_DESC_LEN])
call elpa_generalized_eigenvectors_fc(self, a, b, ev, q, sc_desc, error) call elpa_generalized_eigenvectors_fc(self, a, b, ev, q, sc_desc, .false., error)
end subroutine end subroutine
#endif #endif
......
#if 0 #if 0
subroutine elpa_transform_generalized_& subroutine elpa_transform_generalized_&
&ELPA_IMPL_SUFFIX& &ELPA_IMPL_SUFFIX&
&(self, a, b, sc_desc, error) &(self, a, b, sc_desc, is_already_decomposed, error)
implicit none implicit none
#include "general/precision_kinds.F90" #include "general/precision_kinds.F90"
class(elpa_impl_t) :: self class(elpa_impl_t) :: self
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
MATH_DATATYPE(kind=rck) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols) MATH_DATATYPE(kind=rck) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols)
#endif #endif
integer :: error integer :: error
logical :: is_already_decomposed
integer :: sc_desc(9) integer :: sc_desc(9)
logical, parameter :: do_use_elpa_hermitian_multiply = .true. logical, parameter :: do_use_elpa_hermitian_multiply = .true.
...@@ -19,17 +20,20 @@ ...@@ -19,17 +20,20 @@
MATH_DATATYPE(kind=rck) :: tmp(self%local_nrows, self%local_ncols) MATH_DATATYPE(kind=rck) :: tmp(self%local_nrows, self%local_ncols)
call self%timer_start("transform_generalized()") call self%timer_start("transform_generalized()")
! TODO: why I cannot use self%elpa ??
! B = U^T*U, B<-U if (.not. is_already_decomposed) then
call self%elpa_cholesky_& ! TODO: why I cannot use self%elpa ??
&ELPA_IMPL_SUFFIX& ! B = U^T*U, B<-U
&(b, error) call self%elpa_cholesky_&
if(error .NE. ELPA_OK) return &ELPA_IMPL_SUFFIX&
! B <- inv(U) &(b, error)
call self%elpa_invert_trm_& if(error .NE. ELPA_OK) return
&ELPA_IMPL_SUFFIX& ! B <- inv(U)
&(b, error) call self%elpa_invert_trm_&
if(error .NE. ELPA_OK) return &ELPA_IMPL_SUFFIX&
&(b, error)
if(error .NE. ELPA_OK) return
end if
if(do_use_elpa_hermitian_multiply) then if(do_use_elpa_hermitian_multiply) then
! tmp <- inv(U^T) * A ! tmp <- inv(U^T) * A
......
...@@ -75,6 +75,9 @@ error: define either TEST_ALL_KERNELS or a valid TEST_KERNEL ...@@ -75,6 +75,9 @@ error: define either TEST_ALL_KERNELS or a valid TEST_KERNEL
#endif #endif
#endif #endif
#ifdef TEST_GENERALIZED_DECOMP_EIGENPROBLEM
#define TEST_GENERALIZED_EIGENPROBLEM
#endif
#ifdef TEST_SINGLE #ifdef TEST_SINGLE
# define EV_TYPE real(kind=C_FLOAT) # define EV_TYPE real(kind=C_FLOAT)
...@@ -625,7 +628,17 @@ program test ...@@ -625,7 +628,17 @@ program test
#if defined(TEST_GENERALIZED_EIGENPROBLEM) #if defined(TEST_GENERALIZED_EIGENPROBLEM)
call e%timer_start("e%generalized_eigenvectors()") call e%timer_start("e%generalized_eigenvectors()")
call e%generalized_eigenvectors(a, b, ev, z, sc_desc, error) #if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
call e%timer_start("is_already_decomposed=.false.")
#endif
call e%generalized_eigenvectors(a, b, ev, z, sc_desc, .false., error)
#if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
call e%timer_stop("is_already_decomposed=.false.")
a = as
call e%timer_start("is_already_decomposed=.true.")
call e%generalized_eigenvectors(a, b, ev, z, sc_desc, .true., error)
call e%timer_stop("is_already_decomposed=.true.")
#endif
call e%timer_stop("e%generalized_eigenvectors()") call e%timer_stop("e%generalized_eigenvectors()")
#endif #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