diff --git a/generate_automake_test_programs.py b/generate_automake_test_programs.py
index e288d75f547b689ff4412e41ab4768132eff468c..8558114fa0ffd12e949a88e895339a096e8449ea 100755
--- a/generate_automake_test_programs.py
+++ b/generate_automake_test_programs.py
@@ -45,6 +45,7 @@ test_type_flag = {
"cholesky": "-DTEST_CHOLESKY",
"hermitian_multiply": "-DTEST_HERMITIAN_MULTIPLY",
"generalized" : "-DTEST_GENERALIZED_EIGENPROBLEM",
+ "generalized_decomp": "-DTEST_GENERALIZED_DECOMP_EIGENPROBLEM",
}
layout_flag = {
@@ -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")):
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
if (t == "cholesky" and ((not (m == "toeplitz" or m == "random")) or s == "2stage")):
continue
diff --git a/src/elpa_api.F90 b/src/elpa_api.F90
index 000ba680862bda2c6ba254a4ec99229adb06681c..498f83bbc0f6e814a5286c86c225805ce02fe1a1 100644
--- a/src/elpa_api.F90
+++ b/src/elpa_api.F90
@@ -821,9 +821,10 @@ module elpa_api
!> \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
+ !> \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
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 elpa_constants
import elpa_t
@@ -837,7 +838,7 @@ module elpa_api
#endif
real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
-
+ logical :: is_already_decomposed
integer, optional :: error
end subroutine
end interface
@@ -858,9 +859,10 @@ module elpa_api
!> \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
+ !> \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
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 elpa_constants
import elpa_t
@@ -874,6 +876,7 @@ module elpa_api
#endif
real(kind=c_float) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
+ logical :: is_already_decomposed
integer, optional :: error
end subroutine
@@ -895,9 +898,10 @@ module elpa_api
!> \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
+ !> \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
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 elpa_constants
import elpa_t
@@ -912,6 +916,7 @@ module elpa_api
#endif
real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
+ logical :: is_already_decomposed
integer, optional :: error
end subroutine
@@ -933,9 +938,10 @@ module elpa_api
!> \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
+ !> \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
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 elpa_constants
import elpa_t
@@ -949,6 +955,7 @@ module elpa_api
#endif
real(kind=c_float) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
+ logical :: is_already_decomposed
integer, optional :: error
end subroutine
diff --git a/src/elpa_impl.F90 b/src/elpa_impl.F90
index 37c820cd9453089b691bacbf141e2bd765cf6bdf..24277ecbf0b3257f6f64b376fa37e32c792056a2 100644
--- a/src/elpa_impl.F90
+++ b/src/elpa_impl.F90
@@ -1511,7 +1511,7 @@ module elpa_impl
!> 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
- 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 elpa1_impl
use elpa_utilities, only : error_unit
@@ -1526,13 +1526,14 @@ module elpa_impl
#endif
real(kind=c_double) :: ev(self%na)
integer :: sc_desc(SC_DESC_LEN)
+ logical :: is_already_decomposed
integer, optional :: error
integer :: error_l
integer(kind=c_int) :: solver
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
error = error_l
else if (error_l .ne. ELPA_OK) then
@@ -1585,7 +1586,7 @@ module elpa_impl
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 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
@@ -1614,7 +1615,7 @@ module elpa_impl
!> 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
- 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 elpa1_impl
use elpa_utilities, only : error_unit
@@ -1630,12 +1631,13 @@ module elpa_impl
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
+ logical :: is_already_decomposed
integer :: error_l
integer(kind=c_int) :: solver
#ifdef WANT_SINGLE_PRECISION_REAL
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
error = error_l
else if (error_l .ne. ELPA_OK) then
@@ -1693,7 +1695,7 @@ module elpa_impl
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 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
@@ -1722,7 +1724,7 @@ module elpa_impl
!> 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
- 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 elpa1_impl
use elpa_utilities, only : error_unit
@@ -1739,11 +1741,12 @@ module elpa_impl
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
+ logical :: is_already_decomposed
integer :: error_l
integer(kind=c_int) :: solver
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
error = error_l
else if (error_l .ne. ELPA_OK) then
@@ -1798,7 +1801,7 @@ module elpa_impl
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 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
@@ -1827,7 +1830,7 @@ module elpa_impl
!> 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
- 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 elpa1_impl
use elpa_utilities, only : error_unit
@@ -1844,12 +1847,13 @@ module elpa_impl
integer :: sc_desc(SC_DESC_LEN)
integer, optional :: error
+ logical :: is_already_decomposed
integer :: error_l
integer(kind=c_int) :: solver
#ifdef WANT_SINGLE_PRECISION_COMPLEX
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
error = error_l
else if (error_l .ne. ELPA_OK) then
@@ -1908,7 +1912,7 @@ module elpa_impl
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 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
#endif
diff --git a/src/elpa_impl_template.F90 b/src/elpa_impl_template.F90
index 116a7f02450f8bb9c46a01ce2bcb8530c16a7182..872d327ddaa03739ace70bd426949910a55ed77e 100644
--- a/src/elpa_impl_template.F90
+++ b/src/elpa_impl_template.F90
@@ -1,7 +1,7 @@
#if 0
subroutine elpa_transform_generalized_&
&ELPA_IMPL_SUFFIX&
- &(self, a, b, sc_desc, error)
+ &(self, a, b, sc_desc, is_already_decomposed, error)
implicit none
#include "general/precision_kinds.F90"
class(elpa_impl_t) :: self
@@ -11,6 +11,7 @@
MATH_DATATYPE(kind=rck) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols)
#endif
integer :: error
+ logical :: is_already_decomposed
integer :: sc_desc(9)
logical, parameter :: do_use_elpa_hermitian_multiply = .true.
@@ -19,17 +20,20 @@
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_&
- &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
+
+ if (.not. is_already_decomposed) then
+ ! 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
+ end if
if(do_use_elpa_hermitian_multiply) then
! tmp <- inv(U^T) * A
diff --git a/test/Fortran/test.F90 b/test/Fortran/test.F90
index 026b85d4b983325708de1c1f9d58c6687a574575..42d103a3236721453898668865d93d54f4ff56d4 100644
--- a/test/Fortran/test.F90
+++ b/test/Fortran/test.F90
@@ -75,6 +75,9 @@ error: define either TEST_ALL_KERNELS or a valid TEST_KERNEL
#endif
#endif
+#ifdef TEST_GENERALIZED_DECOMP_EIGENPROBLEM
+#define TEST_GENERALIZED_EIGENPROBLEM
+#endif
#ifdef TEST_SINGLE
# define EV_TYPE real(kind=C_FLOAT)
@@ -625,7 +628,17 @@ program test
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
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()")
#endif