Commit e203c017 authored by Daniel Wortmann's avatar Daniel Wortmann

First attempt to implement a MPI interface to ChASE. Compiles but is untested,...

First attempt to implement a MPI interface to ChASE. Compiles but is untested, hence does not work yet...
parent a8061998
......@@ -4,6 +4,7 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!
! @authors: Miriam Hinzen, Gregor Michalicek
! Added MPI implementation, DW 2018
!--------------------------------------------------------------------------------
MODULE m_chase_diag
#ifdef CPP_CHASE
......@@ -30,6 +31,43 @@ IMPLICIT NONE
end subroutine chase_r
end interface
!MPI
INTERFACE
SUBROUTINE mpi_dchase_init( mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen) BIND( c, name = 'dchase_init' )
USE, INTRINSIC :: iso_c_binding
INTEGER(c_int) :: mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen
END SUBROUTINE mpi_dchase_init
END INTERFACE
INTERFACE
SUBROUTINE mpi_zchase_init( mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen) BIND( c, name = 'zchase_init' )
USE, INTRINSIC :: iso_c_binding
INTEGER(c_int) :: mpi_comm, n, nev, nex, xoff,yoff,xlen,ylen
END SUBROUTINE mpi_zchase_init
END INTERFACE
INTERFACE
SUBROUTINE mpi_chase_r(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'dchase_solve_' )
USE, INTRINSIC :: iso_c_binding
REAL(c_double_complex) :: h(*), v(*)
INTEGER(c_int) :: deg
REAL(c_double) :: ritzv(*), tol
CHARACTER(len=1,kind=c_char) :: mode, opt
END SUBROUTINE mpi_chase_r
END INTERFACE
INTERFACE
SUBROUTINE mpi_chase_c(h, v, ritzv, deg, tol, mode, opt ) BIND( c, name = 'zchase_solve_' )
USE, INTRINSIC :: iso_c_binding
COMPLEX(c_double_complex) :: h(*), v(*)
INTEGER(c_int) :: deg
REAL(c_double) :: ritzv(*), tol
CHARACTER(len=1,kind=c_char) :: mode, opt
END SUBROUTINE mpi_chase_c
END INTERFACE
PRIVATE
INTEGER :: chase_eig_id
......@@ -39,7 +77,7 @@ IMPLICIT NONE
CONTAINS
SUBROUTINE init_chase(mpi,dimension,atoms,kpts,noco,l_real)
USE m_types_mpimat
USE m_types
USE m_types_mpi
USE m_judft
......@@ -65,7 +103,44 @@ IMPLICIT NONE
END IF
END SUBROUTINE init_chase
SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
SUBROUTINE chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
USE m_types_mpimat
USE m_types
USE m_judft
USE iso_c_binding
USE m_eig66_io
!Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
IMPLICIT NONE
CLASS(t_mat), INTENT(INOUT) :: hmat,smat
INTEGER, INTENT(IN) :: ikpt
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: iter
INTEGER, INTENT(INOUT) :: ne
CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: zmat
REAL, INTENT(OUT) :: eig(:)
!Choose serial or parallel solver
SELECT TYPE(hmat)
CLASS is (t_mpimat)
SELECT TYPE(smat)
CLASS is (t_mpimat)
CALL chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
CLASS default
CALL judft_error("Inconsistent matrix setup")
END SELECT
CLASS is (t_mat)
SELECT TYPE(smat)
CLASS is (t_mat)
CALL chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
CLASS default
CALL judft_error("Inconsistent matrix setup")
END SELECT
END SELECT
END SUBROUTINE chase_diag
SUBROUTINE chase_diag_noMPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
USE m_types
USE m_judft
......@@ -211,6 +286,128 @@ IMPLICIT NONE
ENDIF
IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
END SUBROUTINE chase_diag
END SUBROUTINE chase_diag_noMPI
SUBROUTINE chase_diag_MPI(hmat,smat,ikpt,jsp,iter,ne,eig,zmat)
use m_types_mpimat
USE m_types
USE m_judft
USE iso_c_binding
USE m_eig66_io
!Simple driver to solve Generalized Eigenvalue Problem using the ChASE library
IMPLICIT NONE
TYPE(t_mpimat), INTENT(INOUT) :: hmat,smat
INTEGER, INTENT(IN) :: ikpt
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: iter
INTEGER, INTENT(INOUT) :: ne
CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: zmat
REAL, INTENT(OUT) :: eig(:)
INTEGER :: i, j, nev, nex, nbands,xoff,yoff,xlen,ylen,ierr
INTEGER :: info
CLASS(t_mat), ALLOCATABLE :: zMatTemp,chase_mat
REAL, ALLOCATABLE :: eigenvalues(:)
include 'mpif.h'
!Transform to standard problem using SCALAPACK
IF (hmat%l_real) THEN
CALL pdpotrf('U',smat%global_size1,smat%data_r,1,1,smat%blacs_desc,info)
IF (info.NE.0) THEN
WRITE (*,*) 'Error in pdpotrf: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
CALL pdsygst(1,'U',smat%global_size1,hmat%data_r,1,1,hmat%blacs_desc,smat%data_r,1,1,smat%blacs_desc,info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in pdsygst: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
ELSE
CALL pzpotrf('U',smat%global_size1,smat%data_c,1,1,smat%blacs_desc,info)
IF (info.NE.0) THEN
WRITE (*,*) 'Error in pzpotrf: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
CALL pzhegst(1,'U',smat%global_size1,hmat%data_c,1,1,hmat%blacs_desc,smat%data_c,1,1,hmat%blacs_desc,info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in pzhegst: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
END IF
! Now we are ready to set up chase
nev = min(ne,hmat%global_size1)
nex = min(max(nev/4, 45), hmat%global_size1-nev) !dimensioning for workspace
IF (smat%l_real) THEN
CALL mpi_dchase_init( hmat%mpi_com,hmat%global_size1, nev, nex, xoff,yoff,xlen,ylen)
ELSE
CALL mpi_zchase_init( hmat%mpi_com,hmat%global_size1, nev, nex, xoff,yoff,xlen,ylen)
ENDIF
CALL MPI_ALLREDUCE(MPI_IN_PLACE,xlen,1,MPI_INTEGER,MPI_MAX,smat%mpi_com,ierr)
CALL MPI_ALLREDUCE(MPI_IN_PLACE,ylen,1,MPI_INTEGER,MPI_MAX,smat%mpi_com,ierr)
ALLOCATE(t_mpimat::chase_mat)
CALL chase_mat%init(smat%l_real,smat%global_size1,smat%global_size2,smat%mpi_com,.TRUE.,xlen,ylen)
CALL chase_mat%copy(hmat,1,1) !redistribute matrix
ALLOCATE(eigenvalues(nev+nex))
eigenvalues = 0.0
ALLOCATE(t_mpimat::zmatTemp)
CALL zMatTemp%init(hmat%l_real,hmat%global_size1,nev+nex,MPI_COMM_SELF,.TRUE.,hmat%global_size1,nev+nex) !Generate a pseudo-distributed matrix
IF (hmat%l_real) THEN
IF(iter.EQ.1) THEN
CALL chase_r(hmat%data_r, hmat%global_size1, zMatTemp%data_r, eigenvalues, nev, nex, 25, 1e-6, 'R', 'S' )
ELSE
CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
CALL chase_r(hmat%data_r, hmat%global_size1, zMatTemp%data_r, eigenvalues, nev, nex, 25, 1e-6, 'A', 'S' )
END IF
ELSE
IF(iter.EQ.1) THEN
CALL chase_c(hmat%data_c, hmat%global_size1, zMatTemp%data_c, eigenvalues, nev, nex, 25, 1e-6, 'R', 'S' )
ELSE
CALL read_eig(chase_eig_id,ikpt,jsp,neig=nbands,eig=eigenvalues,zmat=zMatTemp)
CALL chase_c(hmat%data_c, hmat%global_size1, zMatTemp%data_c, eigenvalues, nev, nex, 25, 1e-6, 'A', 'S' )
END IF
ENDIF
ne = nev
CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
eigenvalues(:(nev+nex)),zmat=zMatTemp)
CALL zmatTemp%copy(hmat,1,1) !Copy matrix into distributed form
call zmatTemp%free()
! --> recover the generalized eigenvectors z by solving z' = l^t * z
IF (smat%l_real) THEN
CALL pdtrtrs('U','N','N',hmat%global_size1,nev,smat%data_r,1,1,smat%blacs_desc,&
hmat%data_r,1,1,smat%blacs_desc,info)
ELSE
CALL pztrtrs('U','N','N',hmat%global_size1,nev,smat%data_c,1,1,smat%blacs_desc,&
hmat%data_c,1,1,smat%blacs_desc,info)
END IF
IF (info.NE.0) THEN
WRITE (6,*) 'Error in p?trtrs: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
ENDIF
! Redistribute eigvec from ScaLAPACK distribution to each process
! having all eigenvectors corresponding to his eigenvalues as above
!
ALLOCATE(t_mpimat::zmat)
CALL zmat%init(hmat%l_real,hmat%global_size1,hmat%global_size1,hmat%mpi_com,.FALSE.)
CALL zmat%copy(hmat,1,1)
DO i = 1, ne
eig(i) = eigenvalues(i)
END DO
END SUBROUTINE chase_diag_MPI
#endif
END MODULE m_chase_diag
......@@ -58,12 +58,12 @@ MODULE m_types_mat
SUBROUTINE t_mat_init(mat,l_real,matsize1,matsize2,mpi_subcom,l_2d)
SUBROUTINE t_mat_init(mat,l_real,matsize1,matsize2,mpi_subcom,l_2d,nb_x,nb_y)
CLASS(t_mat) :: mat
LOGICAL,INTENT(IN),OPTIONAL:: l_real
INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2
INTEGER,INTENT(IN),OPTIONAL:: mpi_subcom !not needed here, only for allowing overloading this in t_mpimat
LOGICAL,INTENT(IN),OPTIONAL:: l_2d !not needed here either
INTEGER,INTENT(IN),OPTIONAL:: mpi_subcom,nb_x,nb_y !not needed here, only for allowing overloading this in t_mpimat
LOGICAL,INTENT(IN),OPTIONAL:: l_2d !not needed here either
CALL mat%alloc(l_real,matsize1,matsize2)
END SUBROUTINE t_mat_init
......
......@@ -120,31 +120,34 @@ CONTAINS
!! - TRUE: the matrix is a Scalapack BLOCK-CYCLIC distribution
!! - FALSE: the matrix is distributed in a one-dimensional column cyclic distribution with blocksize 1
!! as used in the parallel matrix setup of FLEUR
SUBROUTINE mpimat_init(mat,l_real,matsize1,matsize2,mpi_subcom,l_2d)
SUBROUTINE mpimat_init(mat,l_real,matsize1,matsize2,mpi_subcom,l_2d,nb_x,nb_y)
IMPLICIT NONE
CLASS(t_mpimat) :: mat
INTEGER,INTENT(IN),OPTIONAL :: matsize1,matsize2,mpi_subcom
LOGICAL,INTENT(IN),OPTIONAL :: l_real,l_2d
INTEGER,INTENT(IN),OPTIONAL :: nb_y,nb_x
INTEGER::nb
nb=DEFAULT_BLOCKSIZE
INTEGER::nbx,nby
nbx=DEFAULT_BLOCKSIZE; nby=DEFAULT_BLOCKSIZE
IF (PRESENT(nb_x)) nbx=nb_x
IF (PRESENT(nb_y)) nby=nb_y
IF (.NOT.(PRESENT(matsize1).AND.PRESENT(matsize2).AND.PRESENT(mpi_subcom).AND.PRESENT(l_real).AND.PRESENT(l_2d)))&
CALL judft_error("Optional arguments must be present in mpimat_init")
mat%global_size1=matsize1
mat%global_size2=matsize2
mat%mpi_com=mpi_subcom
CALL priv_create_blacsgrid(mat%mpi_com,l_2d,matsize1,matsize2,nb,&
CALL priv_create_blacsgrid(mat%mpi_com,l_2d,matsize1,matsize2,nbx,nby,&
mat%blacs_ctext,mat%blacs_desc,&
mat%matsize1,mat%matsize2,&
mat%npcol,mat%nprow)
CALL mat%alloc(l_real) !Attention,sizes determined in call to priv_create_blacsgrid
END SUBROUTINE mpimat_init
SUBROUTINE priv_create_blacsgrid(mpi_subcom,l_2d,m1,m2,nb,ictextblacs,sc_desc,local_size1,local_size2,npcol,nprow)
SUBROUTINE priv_create_blacsgrid(mpi_subcom,l_2d,m1,m2,nbc,nbr,ictextblacs,sc_desc,local_size1,local_size2,npcol,nprow)
IMPLICIT NONE
INTEGER,INTENT(IN) :: mpi_subcom
INTEGER,INTENT(IN) :: m1,m2
INTEGER,INTENT(INOUT)::nb
INTEGER,INTENT(INOUT)::nbc,nbr
LOGICAL,INTENT(IN) :: l_2d
INTEGER,INTENT(OUT):: ictextblacs,sc_desc(:)
INTEGER,INTENT(OUT):: local_size1,local_size2
......@@ -179,7 +182,8 @@ CONTAINS
ENDIF
ENDDO distloop
ELSE
nb=1
nbc=1
nbr=MAX(m1,m2)
npcol=np
nprow=1
ENDIF
......@@ -197,9 +201,9 @@ CONTAINS
!
! Now allocate Asca to put the elements of Achi or receivebuffer to
!
myrowssca=(m1-1)/(nb*nprow)*nb+ MIN(MAX(m1-(m1-1)/(nb*nprow)*nb*nprow-nb*myrow,0),nb)
myrowssca=(m1-1)/(nbr*nprow)*nbr+ MIN(MAX(m1-(m1-1)/(nbr*nprow)*nbr*nprow-nbr*myrow,0),nbr)
! Number of rows the local process gets in ScaLAPACK distribution
mycolssca=(m2-1)/(nb*npcol)*nb+ MIN(MAX(m2-(m2-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb)
mycolssca=(m2-1)/(nbc*npcol)*nbc+ MIN(MAX(m2-(m2-1)/(nbc*npcol)*nbc*npcol-nbc*mycol,0),nbc)
!Get BLACS ranks for all MPI ranks
CALL BLACS_PINFO(iamblacs,npblacs) ! iamblacs = local process rank (e.g. myid)
......@@ -237,7 +241,7 @@ CONTAINS
ENDIF
!Create the descriptors
CALL descinit(sc_desc,m1,m2,nb,nb,0,0,ictextblacs,myrowssca,ierr)
CALL descinit(sc_desc,m1,m2,nbr,nbc,0,0,ictextblacs,myrowssca,ierr)
IF (ierr /=0 ) call judft_error('Creation of BLACS descriptor failed')
local_size1=myrowssca
local_size2=mycolssca
......
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