Commit c83ccea4 authored by Gregor Michalicek's avatar Gregor Michalicek

Next step in making ChASE use eigenvectors of the previous iteration

...together with Miriam Hinzen.
parent 6fdd343c
......@@ -32,15 +32,21 @@ IMPLICIT NONE
CONTAINS
SUBROUTINE chase_diag(hmat,smat,ne,eig,zmat)
SUBROUTINE chase_diag(mpi,hmat,smat,ikpt,jsp,chase_eig_id,ne,eig,zmat)
USE m_types
USE m_types_mpi
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_mpi), INTENT(IN) :: mpi
TYPE(t_mat), INTENT(INOUT) :: hmat,smat
INTEGER, INTENT(IN) :: ikpt
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: chase_eig_id
INTEGER, INTENT(INOUT) :: ne
CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: zmat
REAL, INTENT(OUT) :: eig(:)
......@@ -48,14 +54,21 @@ IMPLICIT NONE
INTEGER :: i, j, nev, nex
INTEGER :: info
CLASS(t_Mat), ALLOCATABLE :: zMatTemp
REAL(c_double), ALLOCATABLE :: eigenvalues(:)
REAL(c_double), ALLOCATABLE :: eigenvectors_r(:,:)
COMPLEX(c_double_complex), ALLOCATABLE :: eigenvectors_c(:,:)
ALLOCATE(t_mat::zmat)
CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)
nev = min(ne,hmat%matsize1)
nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
ALLOCATE(eigenvalues(nev+nex))
eigenvalues = 0.0
ALLOCATE(t_mat::zmatTemp)
CALL zMatTemp%alloc(hmat%l_real,hmat%matsize1,nev+nex)
IF (hmat%l_real) THEN
! --> start with Cholesky factorization of b ( so that b = l * l^t)
......@@ -76,14 +89,7 @@ IMPLICIT NONE
! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
nev = min(ne,hmat%matsize1)
nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
ALLOCATE(eigenvectors_r(smat%matsize1,nev+nex))
ALLOCATE(eigenvalues(nev+nex))
eigenvectors_r = 0.0
eigenvalues = 0.0
zMatTemp%data_r = 0.0
do j = 1, hmat%matsize1
do i = 1, j
......@@ -91,15 +97,21 @@ IMPLICIT NONE
end do
end do
! if(first_entry_franza) then
call chase_r(hmat%data_r, hmat%matsize1, eigenvectors_r, eigenvalues, nev, nex, 25, 1e-10, 'R', 'S' )
call chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 25, 1e-10, 'R', 'S' )
! else
! call chase_r(hmat%data_r, hmat%matsize1, eigenvectors_r, eigenvalues, nev, nex, 25, 1e-10, 'A', 'S' )
! load eigenvectors,...
! call chase_r(hmat%data_r, hmat%matsize1, zMatTemp%data_r, eigenvalues, nev, nex, 25, 1e-10, 'A', 'S' )
! end if
ne = nev
CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
eigenvalues(:(nev+nex)),n_start=mpi%n_size,n_end=mpi%n_rank,zmat=zMatTemp)
!TODO: Store eigenvectors array to reuse it in next iteration
! --> recover the generalized eigenvectors z by solving z' = l^t * z
CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,eigenvectors_r,zmat%matsize1,info)
CALL dtrtrs('U','N','N',hmat%matsize1,nev,smat%data_r,smat%matsize1,zMatTemp%data_r,zmat%matsize1,info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in dtrtrs: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
......@@ -107,20 +119,12 @@ IMPLICIT NONE
DO i = 1, ne
DO j = 1, hmat%matsize1
zmat%data_r(j,i) = eigenvectors_r(j,i)
zmat%data_r(j,i) = zMatTemp%data_r(j,i)
END DO
eig(i) = eigenvalues(i)
END DO
!TODO: Store eigenvectors array to reuse it in next iteration
DEALLOCATE(eigenvalues)
DEALLOCATE(eigenvectors_r)
ELSE
! --> start with Cholesky factorization of b ( so that b = l * l^t)
......@@ -141,14 +145,7 @@ IMPLICIT NONE
! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub
nev = min(ne,hmat%matsize1)
nex = min(max(nev/4, 45), hmat%matsize1-nev) !dimensioning for workspace
ALLOCATE(eigenvectors_c(smat%matsize1,nev+nex))
ALLOCATE(eigenvalues(nev+nex))
eigenvectors_c = CMPLX(0.0,0.0)
eigenvalues = 0.0
zMatTemp%data_c = CMPLX(0.0,0.0)
do j = 1, hmat%matsize1
do i = 1, j
......@@ -157,15 +154,18 @@ IMPLICIT NONE
end do
! if(first_entry_franza) then
call chase_c(hmat%data_c, hmat%matsize1, eigenvectors_c, eigenvalues, nev, nex, 25, 1e-10, 'R', 'S' )
call chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 25, 1e-10, 'R', 'S' )
! else
! call chase_c(hmat%data_c, hmat%matsize1, eigenvectors_c, eigenvalues, nev, nex, 25, 1e-10, 'A', 'S' )
! call chase_c(hmat%data_c, hmat%matsize1, zMatTemp%data_c, eigenvalues, nev, nex, 25, 1e-10, 'A', 'S' )
! end if
ne = nev
CALL write_eig(chase_eig_id,ikpt,jsp,nev+nex,nev+nex,&
eigenvalues(:(nev+nex)),n_start=mpi%n_size,n_end=mpi%n_rank,zmat=zMatTemp)
! --> recover the generalized eigenvectors z by solving z' = l^t * z
CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,eigenvectors_c,zmat%matsize1,info)
CALL ztrtrs('U','N','N',hmat%matsize1,nev,smat%data_c,smat%matsize1,zMatTemp%data_c,zmat%matsize1,info)
IF (info.NE.0) THEN
WRITE (6,*) 'Error in ztrtrs: info =',info
CALL juDFT_error("Diagonalization failed",calledby="chase_diag")
......@@ -173,18 +173,11 @@ IMPLICIT NONE
DO i = 1, ne
DO j = 1, hmat%matsize1
zmat%data_c(j,i) = eigenvectors_c(j,i)
zmat%data_c(j,i) = zMatTemp%data_c(j,i)
END DO
eig(i) = eigenvalues(i)
END DO
!TODO: Store eigenvectors array to reuse it in next iteration
DEALLOCATE(eigenvalues)
DEALLOCATE(eigenvectors_c)
ENDIF
IF (info.NE.0) CALL judft_error("Diagonalization via ChASE failed", calledby = 'chase_diag')
END SUBROUTINE chase_diag
......
......@@ -39,22 +39,27 @@ CONTAINS
parallel_solver_available=any((/diag_elpa,diag_elemental,diag_scalapack/)>0)
END FUNCTION parallel_solver_available
SUBROUTINE eigen_diag(hmat,smat,ne,eig,ev)
SUBROUTINE eigen_diag(mpi,hmat,smat,ikpt,jsp,chase_eig_id,ne,eig,ev)
USE m_lapack_diag
USE m_magma
USE m_elpa
USE m_scalapack
USE m_elemental
USE m_chase_diag
USE m_types_mpi
USE m_types_mpimat
IMPLICIT NONE
#ifdef CPP_MPI
include 'mpif.h'
#endif
CLASS(t_mat),INTENT(INOUT) :: smat,hmat
CLASS(t_mat),ALLOCATABLE,INTENT(OUT) :: ev
INTEGER,INTENT(INOUT) :: ne
REAL,INTENT(OUT) :: eig(:)
TYPE(t_mpi), INTENT(IN) :: mpi
CLASS(t_mat), INTENT(INOUT) :: smat,hmat
CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: ev
INTEGER, INTENT(IN) :: ikpt
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: chase_eig_id
INTEGER, INTENT(INOUT) :: ne
REAL, INTENT(OUT) :: eig(:)
!Locals
LOGICAL :: parallel
......@@ -80,7 +85,7 @@ CONTAINS
CALL lapack_diag(hmat,smat,ne,eig,ev)
CASE (diag_chase)
#ifdef CPP_CHASE
CALL chase_diag(hmat,smat,ne,eig,ev)
CALL chase_diag(mpi,hmat,smat,ikpt,jsp,chase_eig_id,ne,eig,ev)
#else
CALL juDFT_error('ChASE eigensolver selected but not available', calledby = 'eigen_diag')
#endif
......
......@@ -18,9 +18,9 @@ CONTAINS
!>
!> The matrices generated and diagonalized here are of type m_mat as defined in m_types_mat.
!>@author D. Wortmann
SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
sym,kpts,DIMENSION, vacuum, input, cell, enpara,banddos, noco, oneD,hybrid,&
it,eig_id,results,inden,v,vx)
SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,sym,kpts,DIMENSION,vacuum,input,&
cell,enpara,banddos,noco,oneD,hybrid,it,eig_id,chase_eig_id,results,inden,v,vx)
USE m_constants, ONLY : pi_const,sfp_const
USE m_types
USE m_apws
......@@ -67,15 +67,16 @@ CONTAINS
#endif
! ..
! .. Scalar Arguments ..
INTEGER,INTENT(IN) :: it
INTEGER,INTENT(INOUT):: eig_id
INTEGER,INTENT(IN) :: it
INTEGER,INTENT(INOUT) :: eig_id
INTEGER,INTENT(INOUT) :: chase_eig_id
! ..
!-odim
!+odim
! ..
! .. Local Scalars ..
INTEGER jsp,nk,nred,ne_all,ne_found
INTEGER ne ,lh0
INTEGER jsp,nk,nred,ne_all,ne_found,nevd,nexd
INTEGER ne,lh0
INTEGER isp,i,j,err
LOGICAL l_wu,l_file,l_real,l_zref
......@@ -121,11 +122,20 @@ CONTAINS
mpi%n_size,layers=vacuum%layers,nstars=vacuum%nstars,ncored=DIMENSION%nstd,&
nsld=atoms%nat,nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,l_mcd=banddos%l_mcd,&
l_orb=banddos%l_orb)
!
IF (juDFT_was_argument("-diag:chase")) THEN
nevd = min(dimension%neigd,dimension%nvd+atoms%nlotot)
nexd = min(max(nevd/4, 45),dimension%nvd+atoms%nlotot-nevd) !dimensioning for workspace
chase_eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,nevd+nexd,kpts%nkpt,DIMENSION%jspd,atoms%lmaxd,&
atoms%nlod,atoms%ntype,atoms%nlotot,noco%l_noco,.TRUE.,l_real,noco%l_soc,.FALSE.,&
mpi%n_size,layers=vacuum%layers,nstars=vacuum%nstars,ncored=DIMENSION%nstd,&
nsld=atoms%nat,nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,l_mcd=banddos%l_mcd,&
l_orb=banddos%l_orb)
END IF
!---> set up and solve the eigenvalue problem
!---> loop over spins
!---> set up k-point independent t(l'm',lm) matrices
!
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,v,mpi,results,DIMENSION,td,ud)
neigBuffer = 0
......@@ -163,7 +173,7 @@ CONTAINS
l_wu=.FALSE.
ne_all=DIMENSION%neigd
if (allocated(zmat)) deallocate(zmat)
CALL eigen_diag(hmat,smat,ne_all,eig,zMat)
CALL eigen_diag(mpi,hmat,smat,nk,jsp,chase_eig_id,ne_all,eig,zMat)
DEALLOCATE(hmat,smat)
!
!---> output results
......
......@@ -101,7 +101,7 @@ CONTAINS
CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
! .. Local Scalars ..
INTEGER:: eig_id, archiveType
INTEGER:: eig_id,chase_eig_id, archiveType
INTEGER:: n,it,ithf
LOGICAL:: l_opti,l_cont,l_qfix, l_wann_inp
REAL :: fermiEnergyTemp, fix
......@@ -249,8 +249,8 @@ CONTAINS
vTemp = vTot
CALL enpara%update(mpi,atoms,vacuum,input,vToT)
CALL eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
sym,kpts,DIMENSION,vacuum,input,cell,enpara,banddos,noco,oneD,hybrid,&
it,eig_id,results,inDen,vTemp,vx)
sym,kpts,DIMENSION,vacuum,input,cell,enpara,banddos,noco,oneD,hybrid,&
it,eig_id,chase_eig_id,results,inDen,vTemp,vx)
vTot%mmpMat = vTemp%mmpMat
!!$ eig_idList(pc) = eig_id
CALL timestop("eigen")
......
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