From 7684a40c7c11ef8aa603179cd71e9cb980250fd6 Mon Sep 17 00:00:00 2001 From: Daniel Wortmann Date: Wed, 4 Jul 2018 13:26:39 +0200 Subject: [PATCH] Added functionality to test diagonalization (started ..) --- diagonalization/eigen_diag.F90 | 76 ++++++++++++++++++++++++++++++++-- eigen/eigen.F90 | 2 +- eigen/eigen_redist_matrix.f90 | 9 +++- global/CMakeLists.txt | 1 + global/matrix_copy.f90 | 58 ++++++++++++++++++++++++++ types/types_mat.F90 | 25 ++++++----- types/types_mpimat.F90 | 15 +++++-- 7 files changed, 166 insertions(+), 20 deletions(-) create mode 100644 global/matrix_copy.f90 diff --git a/diagonalization/eigen_diag.F90 b/diagonalization/eigen_diag.F90 index 2d382a3a..66c485b6 100644 --- a/diagonalization/eigen_diag.F90 +++ b/diagonalization/eigen_diag.F90 @@ -39,7 +39,7 @@ CONTAINS parallel_solver_available=any((/diag_elpa,diag_elemental,diag_scalapack/)>0) END FUNCTION parallel_solver_available - SUBROUTINE eigen_diag(hmat,smat,ikpt,jsp,iter,ne,eig,ev) + SUBROUTINE eigen_diag(mpi,hmat,smat,ikpt,jsp,iter,ne,eig,ev) USE m_lapack_diag USE m_magma USE m_elpa @@ -47,10 +47,12 @@ CONTAINS USE m_elemental USE m_chase_diag USE m_types_mpimat + USE m_matrix_copy IMPLICIT NONE #ifdef CPP_MPI include 'mpif.h' #endif + TYPE(t_mpi), INTENT(IN) :: mpi CLASS(t_mat), INTENT(INOUT) :: smat,hmat CLASS(t_mat), ALLOCATABLE, INTENT(OUT) :: ev INTEGER, INTENT(IN) :: ikpt @@ -59,15 +61,33 @@ CONTAINS INTEGER, INTENT(INOUT) :: ne REAL, INTENT(OUT) :: eig(:) - !Locals + !Locals LOGICAL :: parallel + !For check-mode + TYPE(t_mat) :: s_store,h_store + + SELECT TYPE(smat) CLASS IS (t_mpimat) parallel=.TRUE. CLASS default parallel=.FALSE. END SELECT - + + !Create a copy of the matrix if in test mode + IF (judft_was_argument("-diag:test")) THEN + SELECT TYPE(hmat) + CLASS IS (t_mpimat) + CALL s_store%init(hmat%l_real,hmat%global_size1,hmat%global_size2) + CALL h_store%init(hmat%l_real,hmat%global_size1,hmat%global_size2) + CLASS default + CALL s_store%init(hmat%l_real,hmat%matsize1,hmat%matsize2) + CALL h_store%init(hmat%l_real,hmat%matsize1,hmat%matsize2) + END SELECT + CALL matrix_copy(smat,s_store) + CALL matrix_copy(hmat,h_store) + END IF + CALL timestart("Diagonalization") !Select the solver SELECT CASE (priv_select_solver(parallel)) @@ -90,11 +110,61 @@ CONTAINS CASE (diag_debugout) CALL priv_debug_out(smat,hmat) END SELECT + + !Create test the solution + IF (judft_was_argument("-diag:test")) THEN + CALL priv_test_solution(mpi,ne,s_store,h_store,eig,ev) + CALL judft_error("Diagonalization tested") + END IF CALL timestop("Diagonalization") ! END SUBROUTINE eigen_diag + SUBROUTINE priv_test_solution(mpi,ne,s_store,h_store,eig1,ev) + USE m_types + USE m_lapack_diag + USE m_matrix_copy + IMPLICIT NONE + TYPE(t_mpi),INTENT(in):: mpi + INTEGER,INTENT(INOUT) :: ne + TYPE(t_mat) :: s_store,h_store + REAL,INTENT(in) :: eig1(:) + CLASS(t_mat) :: ev + + + REAL,ALLOCATABLE::eig2(:) + TYPE(t_mat) ::ev1 + CLASS(t_mat),ALLOCATABLE ::ev2 + INTEGER :: i,irank + + + ALLOCATE(eig2(ne)) + CALL ev1%init(s_store%l_real,s_store%matsize1,s_store%matsize2) + CALL matrix_copy(ev,ev1) + + IF (mpi%irank==0) THEN + CALL lapack_diag(h_store,s_store,ne,eig2,ev2) + + OPEN(99,file="diag.compare") + WRITE(99,*) "Eigenvalues" + DO i=1,ne + WRITE(99,*) i,eig1(i),eig2(i) + ENDDO + WRITE(99,*) "Eigenvectors" + DO i=1,ne + IF (ev1%l_real) THEN + WRITE(99,"(i0,20(1x,f10.5))") i,ev1%data_r(1:10,i) + WRITE(99,"(i0,20(1x,f10.5))") i,ev2%data_r(1:10,i) + ELSE + WRITE(99,"(i0,20(1x,f10.5))") i,ev1%data_c(1:10,i) + WRITE(99,"(i0,20(1x,f10.5))") i,ev2%data_c(1:10,i) + END IF + ENDDO + CLOSE(99) + END IF + END SUBROUTINE priv_test_solution + SUBROUTINE priv_debug_out(smat,hmat) USE m_types use m_judft diff --git a/eigen/eigen.F90 b/eigen/eigen.F90 index fa304dfe..eb6010ee 100644 --- a/eigen/eigen.F90 +++ b/eigen/eigen.F90 @@ -183,7 +183,7 @@ CONTAINS end select END IF - CALL eigen_diag(hmat,smat,nk,jsp,iter,ne_all,eig,zMat) + CALL eigen_diag(mpi,hmat,smat,nk,jsp,iter,ne_all,eig,zMat) DEALLOCATE(hmat,smat) ! Output results diff --git a/eigen/eigen_redist_matrix.f90 b/eigen/eigen_redist_matrix.f90 index ec04ee39..67317b0f 100644 --- a/eigen/eigen_redist_matrix.f90 +++ b/eigen/eigen_redist_matrix.f90 @@ -32,10 +32,15 @@ CONTAINS CALL mat_final%init(mat(1,1)%l_real,m,m,mpi%sub_comm,.TRUE.) !here the .true. creates a block-cyclic scalapack distribution !up-up component (or only component in collinear case) + IF (SIZE(mat)==1) THEN + CALL mat_final%move(mat(1,1)) + CALL mat(1,1)%free() + RETURN + ENDIF + CALL mat_final%copy(mat(1,1),1,1) CALL mat(1,1)%free() - IF (SIZE(mat)==1) RETURN - + !down-down component CALL mat_final%copy(mat(2,2),lapw%nv(1)+atoms%nlotot+1,lapw%nv(1)+atoms%nlotot+1) CALL mat(2,2)%free() diff --git a/global/CMakeLists.txt b/global/CMakeLists.txt index a74fb822..eeec1407 100644 --- a/global/CMakeLists.txt +++ b/global/CMakeLists.txt @@ -19,6 +19,7 @@ global/savewigner.f ) set(fleur_F90 ${fleur_F90} global/constants.f90 +global/matrix_copy.F90 global/checkdop.F90 global/checkdopall.f90 global/genMTBasis.f90 diff --git a/global/matrix_copy.f90 b/global/matrix_copy.f90 new file mode 100644 index 00000000..5c6dd8a9 --- /dev/null +++ b/global/matrix_copy.f90 @@ -0,0 +1,58 @@ +!-------------------------------------------------------------------------------- +! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany +! This file is part of FLEUR and available as free software under the conditions +! of the MIT license as expressed in the LICENSE file in more detail. +!-------------------------------------------------------------------------------- +!<@brief + !t_mat_to_packed !> convert to packed-storage matrix PROCEDURE :: clear => t_mat_clear !> set data arrays to zero PROCEDURE :: copy => t_mat_copy !> copy into another t_mat (overloaded for t_mpimat) + PROCEDURE :: move => t_mat_move !> move data into another t_mat (overloaded for t_mpimat) PROCEDURE :: init => t_mat_init !> initalize the matrix(overloaded for t_mpimat) PROCEDURE :: free => t_mat_free !> dealloc the data (overloaded for t_mpimat) PROCEDURE :: add_transpose => t_mat_add_transpose!> add the tranpose/Hermitian conjg. without the diagonal (overloaded for t_mpimat) @@ -219,24 +220,26 @@ MODULE m_types_mat end if end subroutine t_mat_inverse + SUBROUTINE t_mat_move(mat,mat1) + IMPLICIT NONE + CLASS(t_mat),INTENT(INOUT):: mat + CLASS(t_mat),INTENT(INOUT):: mat1 + !Special case, the full matrix is copied. Then use move alloc + IF (mat%l_real) THEN + CALL move_ALLOC(mat1%data_r,mat%data_r) + ELSE + CALL move_ALLOC(mat1%data_c,mat%data_c) + END IF + END SUBROUTINE t_mat_move + SUBROUTINE t_mat_copy(mat,mat1,n1,n2) IMPLICIT NONE CLASS(t_mat),INTENT(INOUT):: mat - CLASS(t_mat),INTENT(INOUT) :: mat1 + CLASS(t_mat),INTENT(IN) :: mat1 INTEGER,INTENT(IN) :: n1,n2 INTEGER:: i1,i2 - IF (n1==1.AND.n2==1.AND.mat%matsize1==mat1%matsize1.AND.mat%matsize2==mat1%matsize2) THEN - !Special case, the full matrix is copied. Then use move alloc - IF (mat%l_real) THEN - CALL move_ALLOC(mat1%data_r,mat%data_r) - ELSE - CALL move_ALLOC(mat1%data_c,mat%data_c) - END IF - RETURN - END IF - i1=mat%matsize1-n1+1 !space available for first dimension i2=mat%matsize2-n1+1 i1=MIN(i1,mat1%matsize1) diff --git a/types/types_mpimat.F90 b/types/types_mpimat.F90 index 742eedc3..a96a9fad 100644 --- a/types/types_mpimat.F90 +++ b/types/types_mpimat.F90 @@ -27,6 +27,7 @@ MODULE m_types_mpimat INTEGER:: npcol,nprow !> the number of columns/rows in the processor grid CONTAINS PROCEDURE,PASS :: copy => mpimat_copy ! mpimat_move ! mpimat_free ! mpimat_init ! mpimat_add_transpose !0) mat%blacs_desc(2)=-1 END IF +#endif END SUBROUTINE mpimat_init SUBROUTINE priv_create_blacsgrid(mpi_subcom,l_2d,m1,m2,nbc,nbr,ictextblacs,sc_desc,local_size1,local_size2,nprow,npcol) -- GitLab