Commit a0ceda61 authored by Pavel Kus's avatar Pavel Kus

C interface for generalized eigenproblem

Added C interface for generalized eigen problem
Extended the functionality of test.c to check generalized eigenproblem,
both with and without already decomposed matrix
parent e93f1c06
......@@ -61,6 +61,34 @@
)(handle, a, ev, q, error)
/*! \brief generic C method for elpa_generalized_eigenvectors
*
* \details
* \param handle handle of the ELPA object, which defines the problem
* \param a float/double float complex/double complex pointer to matrix a
* \param b float/double float complex/double complex pointer to matrix b
* \param ev on return: float/double pointer to eigenvalues
* \param q on return: float/double float complex/double complex pointer to eigenvectors
* \param is_already_decomposed set to 1, if b already decomposed by previous call to elpa_generalized
* \param sc_desc scalapack descriptor
* \param error on return the error code, which can be queried with elpa_strerr()
* \result void
*/
#define elpa_generalized_eigenvectors(handle, a, b, ev, q, sc_desc, is_already_decomposed, error) _Generic((a), \
double*: \
elpa_generalized_eigenvectors_d, \
\
float*: \
elpa_generalized_eigenvectors_f, \
\
double complex*: \
elpa_generalized_eigenvectors_dc, \
\
float complex*: \
elpa_generalized_eigenvectors_fc \
)(handle, a, b, ev, q, sc_desc, is_already_decomposed, error)
/*! \brief generic C method for elpa_eigenvalues
*
* \details
......
......@@ -1569,14 +1569,17 @@ module elpa_impl
endif
end subroutine
!c> void elpa_generalized_eigenvectors_d(elpa_t handle, double *a, double *ev, double *q, int *error);
subroutine elpa_generalized_eigenvectors_d_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, error) &
!c> void elpa_generalized_eigenvectors_d(elpa_t handle, double *a, double *b, double *ev, double *q,
!c> int sc_desc[9], int is_already_decomposed, int *error);
subroutine elpa_generalized_eigenvectors_d_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, is_already_decomposed, error) &
bind(C, name="elpa_generalized_eigenvectors_d")
type(c_ptr), intent(in), value :: handle, a_p, b_p, ev_p, q_p, sc_desc_p
integer(kind=c_int), intent(in), value :: is_already_decomposed
integer(kind=c_int), optional, intent(in) :: error
real(kind=c_double), pointer :: a(:, :), b(:, :), q(:, :), ev(:)
integer(kind=c_int), pointer :: sc_desc(:)
logical :: is_already_decomposed_fortran
type(elpa_impl_t), pointer :: self
call c_f_pointer(handle, self)
......@@ -1585,8 +1588,13 @@ module elpa_impl
call c_f_pointer(ev_p, ev, [self%na])
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])
if(is_already_decomposed .eq. 0) then
is_already_decomposed_fortran = .false.
else
is_already_decomposed_fortran = .true.
end if
call elpa_generalized_eigenvectors_d(self, a, b, ev, q, sc_desc, .false., error)
call elpa_generalized_eigenvectors_d(self, a, b, ev, q, sc_desc, is_already_decomposed_fortran, error)
end subroutine
......@@ -1678,14 +1686,17 @@ module elpa_impl
end subroutine
!c> void elpa_generalized_eigenvectors_f(elpa_t handle, float *a, float *ev, float *q, int *error);
subroutine elpa_generalized_eigenvectors_f_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, error) &
!c> void elpa_generalized_eigenvectors_f(elpa_t handle, float *a, float *b, float *ev, float *q,
!c> int sc_desc[9], int is_already_decomposed, int *error);
subroutine elpa_generalized_eigenvectors_f_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, is_already_decomposed, error) &
bind(C, name="elpa_generalized_eigenvectors_f")
type(c_ptr), intent(in), value :: handle, a_p, b_p, ev_p, q_p, sc_desc_p
integer(kind=c_int), intent(in), value :: is_already_decomposed
integer(kind=c_int), optional, intent(in) :: error
real(kind=c_float), pointer :: a(:, :), b(:, :), q(:, :), ev(:)
integer(kind=c_int), pointer :: sc_desc(:)
logical :: is_already_decomposed_fortran
type(elpa_impl_t), pointer :: self
call c_f_pointer(handle, self)
......@@ -1694,8 +1705,13 @@ module elpa_impl
call c_f_pointer(ev_p, ev, [self%na])
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])
if(is_already_decomposed .eq. 0) then
is_already_decomposed_fortran = .false.
else
is_already_decomposed_fortran = .true.
end if
call elpa_generalized_eigenvectors_f(self, a, b, ev, q, sc_desc, .false., error)
call elpa_generalized_eigenvectors_f(self, a, b, ev, q, sc_desc, is_already_decomposed_fortran, error)
end subroutine
......@@ -1783,15 +1799,18 @@ module elpa_impl
end subroutine
!c> void elpa_generalized_eigenvectors_dc(elpa_t handle, double complex *a, double *ev, double complex *q, int *error);
subroutine elpa_generalized_eigenvectors_dc_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, error) &
!c> void elpa_generalized_eigenvectors_dc(elpa_t handle, double complex *a, double complex *b, double *ev, double complex *q,
!c> int sc_desc[9], int is_already_decomposed, int *error);
subroutine elpa_generalized_eigenvectors_dc_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, is_already_decomposed, error) &
bind(C, name="elpa_generalized_eigenvectors_dc")
type(c_ptr), intent(in), value :: handle, a_p, b_p, ev_p, q_p, sc_desc_p
integer(kind=c_int), intent(in), value :: is_already_decomposed
integer(kind=c_int), optional, intent(in) :: error
complex(kind=c_double_complex), pointer :: a(:, :), b(:, :), q(:, :)
real(kind=c_double), pointer :: ev(:)
integer(kind=c_int), pointer :: sc_desc(:)
logical :: is_already_decomposed_fortran
type(elpa_impl_t), pointer :: self
call c_f_pointer(handle, self)
......@@ -1800,8 +1819,13 @@ module elpa_impl
call c_f_pointer(ev_p, ev, [self%na])
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])
if(is_already_decomposed .eq. 0) then
is_already_decomposed_fortran = .false.
else
is_already_decomposed_fortran = .true.
end if
call elpa_generalized_eigenvectors_dc(self, a, b, ev, q, sc_desc, .false., error)
call elpa_generalized_eigenvectors_dc(self, a, b, ev, q, sc_desc, is_already_decomposed_fortran, error)
end subroutine
......@@ -1894,15 +1918,18 @@ module elpa_impl
end subroutine
!c> void elpa_generalized_eigenvectors_fc(elpa_t handle, float complex *a, float *ev, float complex *q, int *error);
subroutine elpa_generalized_eigenvectors_fc_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, error) &
!c> void elpa_generalized_eigenvectors_fc(elpa_t handle, float complex *a, float complex *b, float *ev, float complex *q,
!c> int sc_desc[9], int is_already_decomposed, int *error);
subroutine elpa_generalized_eigenvectors_fc_c(handle, a_p, b_p, ev_p, q_p, sc_desc_p, is_already_decomposed, error) &
bind(C, name="elpa_generalized_eigenvectors_fc")
type(c_ptr), intent(in), value :: handle, a_p, b_p, ev_p, q_p, sc_desc_p
integer(kind=c_int), intent(in), value :: is_already_decomposed
integer(kind=c_int), optional, intent(in) :: error
complex(kind=c_float_complex), pointer :: a(:, :), b(:, :), q(:, :)
real(kind=c_float), pointer :: ev(:)
integer(kind=c_int), pointer :: sc_desc(:)
logical :: is_already_decomposed_fortran
type(elpa_impl_t), pointer :: self
call c_f_pointer(handle, self)
......@@ -1911,8 +1938,13 @@ module elpa_impl
call c_f_pointer(ev_p, ev, [self%na])
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])
if(is_already_decomposed .eq. 0) then
is_already_decomposed_fortran = .false.
else
is_already_decomposed_fortran = .true.
end if
call elpa_generalized_eigenvectors_fc(self, a, b, ev, q, sc_desc, .false., error)
call elpa_generalized_eigenvectors_fc(self, a, b, ev, q, sc_desc, is_already_decomposed_fortran, error)
end subroutine
#endif
......
......@@ -274,7 +274,8 @@ int main(int argc, char** argv) {
}
/* check the results */
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
if (status !=0){
printf("The computed EVs are not correct !\n");
......
......@@ -278,7 +278,8 @@ int main(int argc, char** argv) {
/* check the results */
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
if (status !=0){
printf("The computed EVs are not correct !\n");
......
......@@ -274,7 +274,8 @@ int main(int argc, char** argv) {
}
/* check the results */
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
if (status !=0){
printf("The computed EVs are not correct !\n");
......
......@@ -277,7 +277,8 @@ int main(int argc, char** argv) {
}
/* check the results */
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
if (status !=0){
printf("The computed EVs are not correct !\n");
......
......@@ -241,9 +241,11 @@ int main(int argc, char** argv) {
/* check the results */
#ifdef DOUBLE_PRECISION_COMPLEX
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#else
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#endif
if (status !=0){
......
......@@ -233,9 +233,12 @@ int main(int argc, char** argv) {
#ifdef DOUBLE_PRECISION_REAL
/* check the results */
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#else
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#endif
if (status !=0){
printf("The computed EVs are not correct !\n");
......
......@@ -236,9 +236,11 @@ int main(int argc, char** argv) {
/* check the results */
#ifdef DOUBLE_PRECISION_COMPLEX
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#else
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#endif
if (status !=0){
printf("The computed EVs are not correct !\n");
......
......@@ -230,9 +230,11 @@ int main(int argc, char** argv) {
/* check the results */
#ifdef DOUBLE_PRECISION_REAL
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#else
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#endif
if (status !=0){
......
......@@ -66,19 +66,39 @@
#error "define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE"
#endif
#ifdef TEST_GENERALIZED_DECOMP_EIGENPROBLEM
#define TEST_GENERALIZED_EIGENPROBLEM
#endif
#ifdef TEST_SINGLE
# define EV_TYPE float
# ifdef TEST_REAL
# define MATRIX_TYPE float
# define PREPARE_MATRIX_RANDOM prepare_matrix_random_real_single_f
# define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_real_single_f
# define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_real_single_f
# define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_real_single_f
# else
# define MATRIX_TYPE complex float
# define PREPARE_MATRIX_RANDOM prepare_matrix_random_complex_single_f
# define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_complex_single_f
# define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_complex_single_f
# define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_complex_single_f
# endif
#else
# define EV_TYPE double
# ifdef TEST_REAL
# define MATRIX_TYPE double
# define PREPARE_MATRIX_RANDOM prepare_matrix_random_real_double_f
# define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_real_double_f
# define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_real_double_f
# define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_real_double_f
# else
# define MATRIX_TYPE complex double
# define PREPARE_MATRIX_RANDOM prepare_matrix_random_complex_double_f
# define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_complex_double_f
# define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_complex_double_f
# define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_complex_double_f
# endif
#endif
......@@ -101,7 +121,7 @@ int main(int argc, char** argv) {
int my_blacs_ctxt, sc_desc[9], info;
/* The Matrix */
MATRIX_TYPE *a, *as, *z;
MATRIX_TYPE *a, *as, *z, *b, *bs;
EV_TYPE *ev;
int error, status;
......@@ -164,18 +184,12 @@ int main(int argc, char** argv) {
as = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
ev = calloc(na, sizeof(EV_TYPE));
#ifdef TEST_REAL
#ifdef TEST_DOUBLE
prepare_matrix_random_real_double_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
#else
prepare_matrix_random_real_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
#endif
#else
#ifdef TEST_DOUBLE
prepare_matrix_random_complex_double_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
#else
prepare_matrix_random_complex_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
#endif
PREPARE_MATRIX_RANDOM(na, myid, na_rows, na_cols, sc_desc, a, z, as);
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
b = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
bs = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
PREPARE_MATRIX_RANDOM_SPD(na, myid, na_rows, na_cols, sc_desc, b, z, bs, nblk, np_rows, np_cols, my_prow, my_pcol);
#endif
if (elpa_init(CURRENT_API_VERSION) != ELPA_OK) {
......@@ -243,27 +257,30 @@ int main(int argc, char** argv) {
if (myid == 0) {
printf("Solver is set to %d \n", value);
}
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
elpa_generalized_eigenvectors(handle, a, b, ev, z, sc_desc, 0, &error);
#if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
//a = as, so that the problem can be solved again
memcpy(a, as, na_rows * na_cols * sizeof(MATRIX_TYPE));
elpa_generalized_eigenvectors(handle, a, b, ev, z, sc_desc, 1, &error);
#endif
#else
/* Solve EV problem */
elpa_eigenvectors(handle, a, ev, z, &error);
#endif
assert_elpa_ok(error);
elpa_deallocate(handle);
elpa_uninit();
/* check the results */
#ifdef TEST_REAL
#ifdef TEST_DOUBLE
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
status = CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol, bs);
#else
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
#endif
#else
#ifdef TEST_DOUBLE
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
#else
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
#endif
status = CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
#endif
if (status !=0){
......@@ -277,6 +294,10 @@ int main(int argc, char** argv) {
free(z);
free(as);
free(ev);
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
free(b);
free(bs);
#endif
#ifdef WITH_MPI
MPI_Finalize();
......
......@@ -234,19 +234,23 @@ int main(int argc, char** argv) {
/* check the results */
#ifdef TEST_REAL
#ifdef TEST_DOUBLE
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
memcpy(a, as, na_rows*na_cols*sizeof(double));
#else
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
memcpy(a, as, na_rows*na_cols*sizeof(float));
#endif
#else
#ifdef TEST_DOUBLE
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
memcpy(a, as, na_rows*na_cols*sizeof(complex double));
#else
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev,
sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
memcpy(a, as, na_rows*na_cols*sizeof(complex float));
#endif
#endif
......
......@@ -197,24 +197,19 @@
#else /* WITH_MPI */
call PRECISION_GEMM(BLAS_TRANS_OR_CONJ,'N',nev,nev,na,ONE,z,na,z,na,ZERO,tmp1,na)
#endif /* WITH_MPI */
!TODO for the C interface, not all information is passed (zeros instead)
!TODO than this part of the test cannot be done
!TODO either we will not have this part of test at all, or it will be improved
if(nblk > 0) then
! First check, whether the elements on diagonal are 1 .. "normality" of the vectors
err = 0.0_rk
do i=1, nev
if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
err = max(err, abs(tmp1(rowLocal,colLocal) - 1.0_rk))
endif
end do
! First check, whether the elements on diagonal are 1 .. "normality" of the vectors
err = 0.0_rk
do i=1, nev
if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
err = max(err, abs(tmp1(rowLocal,colLocal) - 1.0_rk))
endif
end do
#ifdef WITH_MPI
call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else /* WITH_MPI */
errmax = err
errmax = err
#endif /* WITH_MPI */
if (myid==0) print *,'Maximal error in eigenvector lengths:',errmax
end if
if (myid==0) print *,'Maximal error in eigenvector lengths:',errmax
! Second, find the maximal error in the whole Z**T * Z matrix (its diference from identity matrix)
! Initialize tmp2 to unit matrix
......@@ -250,28 +245,26 @@
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
!c> int check_correctness_evp_numeric_residuals_real_double_f(int na, int nev, int na_rows, int na_cols,
!c> double *as, double *z, double *ev,
!c> int sc_desc[9], int myid);
!c> double *as, double *z, double *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol);
#else
!c> int check_correctness_evp_numeric_residuals_real_single_f(int na, int nev, int na_rows, int na_cols,
!c> float *as, float *z, float *ev,
!c> int sc_desc[9], int myid);
!c> float *as, float *z, float *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol);
#endif
#endif /* REALCASE */
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
!c> int check_correctness_evp_numeric_residuals_complex_double_f(int na, int nev, int na_rows, int na_cols,
!c> complex double *as, complex double *z, double *ev,
!c> int sc_desc[9], int myid);
!c> complex double *as, complex double *z, double *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol);
#else
!c> int check_correctness_evp_numeric_residuals_complex_single_f(int na, int nev, int na_rows, int na_cols,
!c> complex float *as, complex float *z, float *ev,
!c> int sc_desc[9], int myid);
!c> complex float *as, complex float *z, float *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol);
#endif
#endif /* COMPLEXCASE */
......@@ -279,7 +272,7 @@ function check_correctness_evp_numeric_residuals_&
&MATH_DATATYPE&
&_&
&PRECISION&
&_f (na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid) result(status) &
&_f (na, nev, na_rows, na_cols, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol) result(status) &
bind(C,name="check_correctness_evp_numeric_residuals_&
&MATH_DATATYPE&
&_&
......@@ -292,22 +285,82 @@ function check_correctness_evp_numeric_residuals_&
#include "../../src/general/precision_kinds.F90"
integer(kind=c_int) :: status
integer(kind=c_int), value :: na, nev, myid, na_rows, na_cols
integer(kind=c_int), value :: na, nev, myid, na_rows, na_cols, nblk, np_rows, np_cols, my_prow, my_pcol
MATH_DATATYPE(kind=rck) :: as(1:na_rows,1:na_cols), z(1:na_rows,1:na_cols)
real(kind=rck) :: ev(1:na)
integer(kind=c_int) :: sc_desc(1:9)
! TODO: I did not want to add all the variables to the C interface as well
! TODO: I think that we should find a better way to pass this information
! TODO: to all the functions anyway (get it from sc_desc, pass elpa_t, etc..)
status = check_correctness_evp_numeric_residuals_&
&MATH_DATATYPE&
&_&
&PRECISION&
& (na, nev, as, z, ev, sc_desc, 0, myid, 0, 0, 0, 0)
& (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
end function
!---- variant for the generalized eigenproblem
!---- unlike in Fortran, we cannot use optional parameter
!---- we thus define a different function
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
!c> int check_correctness_evp_gen_numeric_residuals_real_double_f(int na, int nev, int na_rows, int na_cols,
!c> double *as, double *z, double *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol,
!c> double *bs);
#else
!c> int check_correctness_evp_gen_numeric_residuals_real_single_f(int na, int nev, int na_rows, int na_cols,
!c> float *as, float *z, float *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol,
!c> float *bs);
#endif
#endif /* REALCASE */
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
!c> int check_correctness_evp_gen_numeric_residuals_complex_double_f(int na, int nev, int na_rows, int na_cols,
!c> complex double *as, complex double *z, double *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol,
!c> complex double *bs);
#else
!c> int check_correctness_evp_gen_numeric_residuals_complex_single_f(int na, int nev, int na_rows, int na_cols,
!c> complex float *as, complex float *z, float *ev, int sc_desc[9],
!c> int nblk, int myid, int np_rows, int np_cols, int my_prow, int my_pcol,
!c> complex float *bs);
#endif
#endif /* COMPLEXCASE */
function check_correctness_evp_gen_numeric_residuals_&
&MATH_DATATYPE&
&_&
&PRECISION&
&_f (na, nev, na_rows, na_cols, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol, bs) result(status) &
bind(C,name="check_correctness_evp_gen_numeric_residuals_&
&MATH_DATATYPE&
&_&
&PRECISION&
&_f")
use iso_c_binding
implicit none