Added functionality to test diagonalization (started ..)

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

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

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()

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

global/matrix_copy.f90

+!--------------------------------------------------------------------------------
+! 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
+  !<This matrix_copy module is needed to copy a distributed into a non-distributed matrix or vice versa
+  !<It calls the usual matrix-copy routines in other cases
+MODULE m_matrix_copy
+  IMPLICIT NONE
+CONTAINS
+  SUBROUTINE matrix_copy(mat_in,mat_out)
+    USE m_types
+    USE m_types_mpimat
+
+    CLASS(t_mat),INTENT(IN)   ::mat_in
+    CLASS(t_mat),INTENT(INOUT)::mat_out
+
+    TYPE(t_mpimat)::tmp_mat
+
+#ifdef CPP_MPI
+    INCLUDE 'mpif.h'
+#else
+    INTEGER:: MPI_COMM_SELF
+#endif
+    
+    SELECT TYPE(mat_in)
+    CLASS is(t_mpimat)
+       SELECT TYPE(mat_out)
+       CLASS is(t_mpimat)
+          CALL mat_out%copy(mat_in,1,1)
+       CLASS is(t_mat)
+          !Copy from t_mpimat to t_mat
+          CALL tmp_mat%init(mat_in%l_real,mat_in%global_size1,mat_in%global_size2,MPI_COMM_SELF)
+          CALL tmp_mat%copy(mat_in,1,1) !redistribute to single matrix
+          IF (tmp_mat%l_real) THEN
+             CALL move_ALLOC(tmp_mat%data_r,mat_out%data_r)
+          ELSE
+             CALL move_ALLOC(tmp_mat%data_r,mat_out%data_r)
+          ENDIF
+       END SELECT
+    CLASS is (t_mat)
+       SELECT TYPE(mat_out)
+       CLASS is(t_mpimat)
+          !Copy from non-distributed t_mat to t_mpimat
+          CALL tmp_mat%init(mat_in%l_real,mat_in%matsize1,mat_in%matsize1,MPI_COMM_SELF)
+          IF (tmp_mat%l_real) THEN
+             tmp_mat%data_r=mat_in%data_r
+          ELSE
+             tmp_mat%data_c=mat_in%data_c
+          ENDIF
+          CALL mat_out%copy(tmp_mat,1,1)
+       CLASS is(t_mat)
+          CALL mat_out%copy(mat_in,1,1)
+       END SELECT
+    END SELECT
+  END SUBROUTINE matrix_copy
+END MODULE m_matrix_copy

types/types_mat.F90

@@ -23,6 +23,7 @@ MODULE m_types_mat
      PROCEDURE        :: to_packed=>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)

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     !<overwriten from t_mat, also performs redistribution
+     PROCEDURE,PASS   :: move => mpimat_move     !<overwriten from t_mat, also performs redistribution
      PROCEDURE,PASS   :: free => mpimat_free     !<overwriten from t_mat, takes care of blacs-grids
      PROCEDURE,PASS   :: init => mpimat_init     !<overwriten from t_mat, also calls alloc in t_mat
      PROCEDURE,PASS   :: add_transpose => mpimat_add_transpose !<overwriten from t_mat
@@ -178,7 +179,7 @@ CONTAINS
   SUBROUTINE mpimat_copy(mat,mat1,n1,n2)
     IMPLICIT NONE
     CLASS(t_mpimat),INTENT(INOUT)::mat
-    CLASS(t_mat),INTENT(INOUT)   ::mat1
+    CLASS(t_mat),INTENT(IN)      ::mat1
     INTEGER,INTENT(IN) ::n1,n2
 #ifdef CPP_SCALAPACK
     SELECT TYPE(mat1)
@@ -193,7 +194,14 @@ CONTAINS
     END SELECT
 #endif    
   END SUBROUTINE mpimat_copy
-  
+
+  SUBROUTINE mpimat_move(mat,mat1)
+    IMPLICIT NONE
+    CLASS(t_mpimat),INTENT(INOUT)::mat
+    CLASS(t_mat),INTENT(INOUT)   ::mat1
+    CALL mat%copy(mat1,1,1)
+  END SUBROUTINE mpimat_move
+
   SUBROUTINE mpimat_free(mat)
     IMPLICIT NONE
     CLASS(t_mpimat),INTENT(INOUT) :: mat
@@ -217,7 +225,7 @@ CONTAINS
     INTEGER,INTENT(IN),OPTIONAL :: matsize1,matsize2,mpi_subcom
     LOGICAL,INTENT(IN),OPTIONAL :: l_real,l_2d
     INTEGER,INTENT(IN),OPTIONAL :: nb_y,nb_x
-    
+#ifdef CPP_SCALAPACK    
     INTEGER::nbx,nby,irank,ierr
     include 'mpif.h'
     nbx=DEFAULT_BLOCKSIZE; nby=DEFAULT_BLOCKSIZE
@@ -238,6 +246,7 @@ CONTAINS
        CALL MPI_COMM_RANK(mpi_subcom,irank,ierr)
        IF (irank>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)