Commit c9f2e983 authored by Daniel Wortmann's avatar Daniel Wortmann

Fixed ELPA for version 2016.05.003

parent 0ee8b28e
......@@ -6,5 +6,5 @@ LINK_LIBRARIES ${FLEUR_LIBRARIES}
message("ELPA Library found:${FLEUR_USE_ELPA}")
if (FLEUR_USE_ELPA)
set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_ELPA")
endif()
\ No newline at end of file
set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_ELPA" "CPP_ELPA_201605003")
endif()
......@@ -163,9 +163,9 @@ CONTAINS
#ifdef CPP_ELPA
CASE (diag_elpa)
IF (hamovlp%l_real) THEN
CALL elpa(lapw%nmat,n,SUB_COMM,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,eig,ne_found)
CALL elpa_diag(lapw%nmat,n,SUB_COMM,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,eig,ne_found)
ELSE
CALL elpa(lapw%nmat,n,SUB_COMM,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c,eig,ne_found)
CALL elpa_diag(lapw%nmat,n,SUB_COMM,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c,eig,ne_found)
ENDIF
#endif
#ifdef CPP_ELEMENTAL
......
......@@ -5,20 +5,21 @@
!--------------------------------------------------------------------------------
MODULE m_elpa
PRIVATE
interface elpa
#ifdef CPP_ELPA
interface elpa_diag
module procedure elpa_r,elpa_c
end interface elpa
end interface elpa_diag
!Module to call elpa library for parallel diagonalization
!uses ssubredist1/2 for redistribution
PUBLIC elpa
PUBLIC elpa_diag
CONTAINS
! First the real version of the code
#define CPP_CHOLESKY cholesky_real
#define CPP_invert_trm invert_trm_real
#define CPP_solve_evp_1stage solve_evp_real_1stage
#define CPP_solve_evp solve_evp_real
#define CPP_solve_evp_2stage solve_evp_real_2stage
#define CPP_REALCOMPLEX real
......@@ -55,6 +56,7 @@ CONTAINS
! Now the complex version
#define CPP_CHOLESKY cholesky_complex
#define CPP_invert_trm invert_trm_complex
#define CPP_solve_evp_1stage solve_evp_complex_1stage
#define CPP_solve_evp solve_evp_complex
#define CPP_solve_evp_2stage solve_evp_complex_2stage
#define CPP_REALCOMPLEX complex
......@@ -190,7 +192,7 @@ CONTAINS
!Create communicators for ELPA
!print *,"creating ELPA comms"
#ifdef CPP_ELPA_NEW
#if (defined CPP_ELPA_201605003)||defined(CPP_ELPA_NEW)
ierr=get_elpa_row_col_comms(mpi_subcom, myrowblacs, mycolblacs,mpi_comm_rows, mpi_comm_cols)
#else
CALL get_elpa_row_col_comms(mpi_subcom, myrowblacs, mycolblacs,mpi_comm_rows, mpi_comm_cols)
......@@ -202,7 +204,8 @@ CONTAINS
IF (ierr /=0 ) CALL juDFT_error('descinit1 failed',calledby='elpa')
END SUBROUTINE priv_create_blacsgrid
#else
LOGICAL :: elpa_used=.false.
#endif
end MODULE m_elpa
......@@ -4,9 +4,9 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
USE m_juDFT
USE elpa1
USE m_subredist1
USE m_subredist2
USE elpa1
#ifdef CPP_ELPA2
USE elpa2
#endif
......@@ -84,7 +84,10 @@
! Please note: cholesky_complex/invert_trm_complex are not trimmed for speed.
! The only reason having them is that the Scalapack counterpart
! PDPOTRF very often fails on higher processor numbers for unknown reasons!
#ifdef CPP_ELPA_NEW
#ifdef CPP_ELPA_201605003
ok=CPP_CHOLESKY (m,bsca,SIZE(bsca,1),nb,mycolssca,mpi_comm_rows,mpi_comm_cols,.false.)
ok=CPP_invert_trm(m,bsca,SIZE(bsca,1),nb,mycolssca,mpi_comm_rows,mpi_comm_cols,.false.)
#elif defined CPP_ELPA_NEW
CALL CPP_CHOLESKY (m,bsca,SIZE(bsca,1),nb,mycolssca,mpi_comm_rows,mpi_comm_cols,.false.,ok)
CALL CPP_invert_trm(m,bsca,SIZE(bsca,1),nb,mycolssca,mpi_comm_rows,mpi_comm_cols,.false.,ok)
#else
......@@ -106,15 +109,23 @@
n_row = numroc (n_col, nb, myrow, 0, nprow)
asca(n_row+1:myrowssca,i) = eigvec(n_row+1:myrowssca,i)
ENDDO
#ifdef CPP_ELPA_201605003
ok=CPP_mult ('U', 'L',m, m,bsca,myrowssca,asca,SIZE(asca,1),nb, mpi_comm_rows, mpi_comm_cols,eigvec,myrowssca)
#else
CALL CPP_mult ('U', 'L',m, m,bsca,myrowssca,asca,SIZE(asca,1),nb, mpi_comm_rows, mpi_comm_cols,eigvec,myrowssca)
#endif
! 2b. tmp2 = eigvec**T
CALL CPP_transpose(m,m,CPP_ONE,eigvec,1,1,sc_desc,CPP_ZERO,tmp2,1,1,sc_desc)
! 2c. A = U**-T * tmp2 ( = U**-T * Aorig * U**-1 )
#ifdef CPP_ELPA_201605003
ok=CPP_mult ('U', 'U', m, m, bsca, SIZE(bsca,1), tmp2,&
SIZE(tmp2,1),nb, mpi_comm_rows, mpi_comm_cols, asca, SIZE(asca,1))
#else
CALL CPP_mult ('U', 'U', m, m, bsca, SIZE(bsca,1), tmp2,&
SIZE(tmp2,1),nb, mpi_comm_rows, mpi_comm_cols, asca, SIZE(asca,1))
#endif
! A is only set in the upper half, solve_evp_real needs a full matrix
! Set lower half from upper half
......@@ -133,7 +144,15 @@
! 3. Calculate eigenvalues/eigenvectors of U**-T * A * U**-1
! Eigenvectors go to eigvec
num2=num
#ifdef CPP_ELPA_NEW
#ifdef CPP_ELPA_201605003
#ifdef CPP_ELPA2
ok=CPP_solve_evp_2stage(m,num2,asca,SIZE(asca,1),&
eig2,eigvec,SIZE(asca,1), nb,mycolssca, mpi_comm_rows, mpi_comm_cols,sub_comm)
#else
ok=CPP_solve_evp_1stage(m, num2,asca,SIZE(asca,1),&
eig2,eigvec, SIZE(asca,1), nb,mycolssca,mpi_comm_rows, mpi_comm_cols)
#endif
#elif defined CPP_ELPA_NEW
#ifdef CPP_ELPA2
err=CPP_solve_evp_2stage(m,num2,asca,SIZE(asca,1),&
eig2,eigvec,SIZE(asca,1), nb,mycolssca, mpi_comm_rows, mpi_comm_cols,sub_comm)
......@@ -156,9 +175,13 @@
! mult_ah_b_complex needs the transpose of U**-1, thus tmp2 = (U**-1)**T
CALL CPP_transpose(m,m,CPP_ONE,bsca,1,1,sc_desc,CPP_ZERO,tmp2,1,1,sc_desc)
#ifdef CPP_ELPA_201605003
ok= CPP_mult ('L', 'N',m, num2, tmp2, SIZE(asca,1),&
eigvec, SIZE(asca,1),nb,mpi_comm_rows, mpi_comm_cols, asca, SIZE(asca,1))
#else
CALL CPP_mult ('L', 'N',m, num2, tmp2, SIZE(asca,1),&
eigvec, SIZE(asca,1),nb,mpi_comm_rows, mpi_comm_cols, asca, SIZE(asca,1))
#endif
! END of ELPA stuff
CALL BLACS_GRIDEXIT(ictextblacs,ierr)
CALL MPI_COMM_FREE(mpi_comm_rows,ierr)
......
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