diff --git a/generate_automake_test_programs.py b/generate_automake_test_programs.py index 4dfc0e8a27b0df803d2a1d9c48bbf39218a7262f..564eb5605d8065d9b176df28b96b9edf47b405af 100755 --- a/generate_automake_test_programs.py +++ b/generate_automake_test_programs.py @@ -133,6 +133,10 @@ for m, g, t, p, d, s, l in product( if (t == "qr" and (s == "1stage" or d == "complex")): continue + #TODO: this does not work at the moment + if(t == "generalized" and (l == "all_layouts" or s == "2stage")): + continue + create_test(m, g, t, p, d, s, l, "fortran") diff --git a/src/elpa_impl.F90 b/src/elpa_impl.F90 index 6b3160cac931dd45622d503a2243d5f554cd35a8..c25986952e24c4f9fbe89ee9b783cae6e6f95871 100644 --- a/src/elpa_impl.F90 +++ b/src/elpa_impl.F90 @@ -125,12 +125,16 @@ module elpa_impl procedure, public :: associate_int => elpa_associate_int !< public method to set some pointers procedure, private :: elpa_transform_generalized_d + procedure, private :: elpa_transform_back_generalized_d procedure, private :: elpa_transform_generalized_dc + procedure, private :: elpa_transform_back_generalized_dc #ifdef WANT_SINGLE_PRECISION_REAL procedure, private :: elpa_transform_generalized_f + procedure, private :: elpa_transform_back_generalized_f #endif #ifdef WANT_SINGLE_PRECISION_COMPLEX procedure, private :: elpa_transform_generalized_fc + procedure, private :: elpa_transform_back_generalized_fc #endif end type elpa_impl_t @@ -1236,6 +1240,11 @@ module elpa_impl logical :: success_l call self%elpa_transform_generalized_d(a, b, sc_desc, error_l) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!" + endif call self%get("solver", solver) if (solver .eq. ELPA_SOLVER_1STAGE) then @@ -1257,6 +1266,13 @@ module elpa_impl else if (.not. success_l) then write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!" endif + + call self%elpa_transform_back_generalized_d(b, q, sc_desc, error_l) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!" + endif end subroutine !c> void elpa_generalized_eigenvectors_d(elpa_t handle, double *a, double *ev, double *q, int *error); @@ -1321,10 +1337,18 @@ module elpa_impl integer :: sc_desc(SC_DESC_LEN) integer, optional :: error + 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) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!" + endif + call self%get("solver",solver) if (solver .eq. ELPA_SOLVER_1STAGE) then success_l = elpa_solve_evp_real_1stage_single_impl(self, a, ev, q) @@ -1345,6 +1369,13 @@ module elpa_impl else if (.not. success_l) then write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!" endif + + call self%elpa_transform_back_generalized_f(b, q, sc_desc, error_l) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!" + endif #else print *,"This installation of the ELPA library has not been build with single-precision support" error = ELPA_ERROR @@ -1415,9 +1446,17 @@ module elpa_impl integer :: sc_desc(SC_DESC_LEN) integer, optional :: error + integer :: error_l integer(kind=c_int) :: solver logical :: success_l + call self%elpa_transform_generalized_dc(a, b, sc_desc, error_l) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!" + endif + call self%get("solver", solver) if (solver .eq. ELPA_SOLVER_1STAGE) then success_l = elpa_solve_evp_complex_1stage_double_impl(self, a, ev, q) @@ -1438,6 +1477,13 @@ module elpa_impl else if (.not. success_l) then write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!" endif + + call self%elpa_transform_back_generalized_dc(b, q, sc_desc, error_l) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!" + endif end subroutine @@ -1505,10 +1551,18 @@ module elpa_impl integer :: sc_desc(SC_DESC_LEN) integer, optional :: error + 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) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_generalized() and you did not check for errors!" + endif + call self%get("solver", solver) if (solver .eq. ELPA_SOLVER_1STAGE) then success_l = elpa_solve_evp_complex_1stage_single_impl(self, a, ev, q) @@ -1529,6 +1583,13 @@ module elpa_impl else if (.not. success_l) then write(error_unit,'(a)') "ELPA: Error in solve() and you did not check for errors!" endif + + call self%elpa_transform_back_generalized_fc(b, q, sc_desc, error_l) + if (present(error)) then + error = error_l + else if (error_l .ne. ELPA_OK) then + write(error_unit,'(a)') "ELPA: Error in transform_back_generalized() and you did not check for errors!" + endif #else print *,"This installation of the ELPA library has not been build with single-precision support" error = ELPA_ERROR diff --git a/src/elpa_impl_template.F90 b/src/elpa_impl_template.F90 index 80172e6353960959b56033e6f814b9cfcfde5ea8..f7ddd76c180a68e5a39e47b78ee976502e27d12e 100644 --- a/src/elpa_impl_template.F90 +++ b/src/elpa_impl_template.F90 @@ -55,3 +55,36 @@ end subroutine + + subroutine elpa_transform_back_generalized_& + &ELPA_IMPL_SUFFIX& + &(self, b, q, sc_desc, error) + implicit none +#include "general/precision_kinds.F90" + class(elpa_impl_t) :: self +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=rck) :: b(self%local_nrows, *), q(self%local_nrows, *) +#else + MATH_DATATYPE(kind=rck) :: b(self%local_nrows, self%local_ncols), q(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: part of eigenvectors only +#ifdef WITH_MPI + ! Q <= inv(U) * Q + call p& + &BLAS_CHAR& + &trmm("L", "U", "N", "N", self%na, self%na, & + ONE, b, 1, 1, sc_desc, q, 1, 1, sc_desc) +#else + call BLAS_CHAR& + &trmm("L", "U", "N", "N", self%na, self%na, & + ONE, b, self%na, q, self%na) +#endif + + end subroutine + diff --git a/test/Fortran/test.F90 b/test/Fortran/test.F90 index 259e4474291710f24fc68f483c3f4907fd57a7ad..5ce5e2be7ba3afebb164cd743a820bfa12e341e6 100644 --- a/test/Fortran/test.F90 +++ b/test/Fortran/test.F90 @@ -344,8 +344,8 @@ program test !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) + diagonalElement = 2.546_rk * ONE + subdiagonalElement = ZERO call prepare_matrix_toeplitz(na, diagonalElement, subdiagonalElement, & d, sd, ds, sds, b, bs, nblk, np_rows, & np_cols, my_prow, my_pcol) @@ -518,7 +518,12 @@ program test ! status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol) !#elif defined(TEST_MATRIX_RANDOM) if (nev .ge. 1) then +#if defined(TEST_GENERALIZED_EIGENPROBLEM) + status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, & + my_pcol, bs) +#else status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol) +#endif else ! zero eigenvectors and no analytic test => toeplitz status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, & diff --git a/test/shared/test_check_correctness_template.F90 b/test/shared/test_check_correctness_template.F90 index 29f11152f03ac88fc4311a7df5b73f16663c9699..bb429f8492c1adc78d73b930771693d10b898ac8 100644 --- a/test/shared/test_check_correctness_template.F90 +++ b/test/shared/test_check_correctness_template.F90 @@ -45,12 +45,13 @@ &MATH_DATATYPE& &_& &PRECISION& - & (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol) result(status) + & (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol, bs) result(status) implicit none #include "../../src/general/precision_kinds.F90" integer(kind=ik) :: status integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol - MATH_DATATYPE(kind=rck), intent(in) :: as(:,:), z(:,:) + MATH_DATATYPE(kind=rck), intent(in) :: as(:,:), z(:,:) + MATH_DATATYPE(kind=rck), intent(in), optional :: bs(:,:) real(kind=rk) :: ev(:) MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2 MATH_DATATYPE(kind=rck) :: xc @@ -76,10 +77,13 @@ real(kind=rk), parameter :: tol_res_real_single = 3e-3_rk real(kind=rk), parameter :: tol_res_complex_double = 5e-12_rk real(kind=rk), parameter :: tol_res_complex_single = 3e-3_rk - real(kind=rk), parameter :: tol_res = tol_res_& + real(kind=rk) :: tol_res = tol_res_& &MATH_DATATYPE& &_& &PRECISION + ! precision of generalized problem is lower + real(kind=rk), parameter :: generalized_penalty = 10.0_rk + ! tolerance for the orthogonality test for different math type/precision setups real(kind=rk), parameter :: tol_orth_real_double = 5e-12_rk real(kind=rk), parameter :: tol_orth_real_single = 9e-4_rk @@ -90,34 +94,50 @@ &_& &PRECISION + if(present(bs)) then + tol_res = generalized_penalty * tol_res + endif status = 0 ! 1. Residual (maximum of || A*Zi - Zi*EVi ||) - ! tmp1 = A * Z - ! as is original stored matrix, Z are the EVs -#ifdef WITH_MPI - call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, as, 1, 1, sc_desc, & - z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc) -#else /* WITH_MPI */ - call PRECISION_GEMM('N','N',na,nev,na,ONE,as,na,z,na,ZERO,tmp1,na) -#endif /* WITH_MPI */ - - - ! tmp2 = Zi*EVi - tmp2(:,:) = z(:,:) + ! tmp1 = Zi*EVi + tmp1(:,:) = z(:,:) do i=1,nev xc = ev(i) #ifdef WITH_MPI call p& &BLAS_CHAR& - &scal(na, xc, tmp2, 1, i, sc_desc, 1) + &scal(na, xc, tmp1, 1, i, sc_desc, 1) #else /* WITH_MPI */ call BLAS_CHAR& - &scal(na,xc,tmp2(:,i),1) + &scal(na,xc,tmp1(:,i),1) #endif /* WITH_MPI */ enddo + ! for generalized EV problem, multiply by bs as well + ! tmp2 = B * tmp1 + if(present(bs)) then +#ifdef WITH_MPI + call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, bs, 1, 1, sc_desc, & + tmp1, 1, 1, sc_desc, ZERO, tmp2, 1, 1, sc_desc) +#else /* WITH_MPI */ + call PRECISION_GEMM('N','N',na,nev,na,ONE,bs,na,tmp1,na,ZERO,tmp2,na) +#endif /* WITH_MPI */ + else + ! normal eigenvalue problem .. no need to multiply + tmp2(:,:) = tmp1(:,:) + end if + + ! tmp1 = A * Z + ! as is original stored matrix, Z are the EVs +#ifdef WITH_MPI + call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, as, 1, 1, sc_desc, & + z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc) +#else /* WITH_MPI */ + call PRECISION_GEMM('N','N',na,nev,na,ONE,as,na,z,na,ZERO,tmp1,na) +#endif /* WITH_MPI */ + ! tmp1 = A*Zi - Zi*EVi tmp1(:,:) = tmp1(:,:) - tmp2(:,:) @@ -166,7 +186,9 @@ endif ! 2. Eigenvector orthogonality - + !TODO for the generalized eigenvector problem, the orthogonality test has to be altered + !TODO at the moment, it is skipped + if(.not. present(bs)) then ! tmp1 = Z**T * Z tmp1 = 0 #ifdef WITH_MPI @@ -223,6 +245,8 @@ status = 1 endif endif + + endif ! skiping test of orthogonality for generalized eigenproblem end function