Modified handling of matrices. Try to ensure that blacs contexts are properly freed.

eigen/eigen.F90

@@ -190,6 +190,7 @@ CONTAINS
           END IF
 
           CALL eigen_diag(mpi,hmat,smat,nk,jsp,iter,ne_all,eig,zMat)
+          CALL smat%free()
           DEALLOCATE(hmat,smat, stat=dealloc_stat, errmsg=errmsg)
           if(dealloc_stat /= 0) call juDFT_error("deallocate failed for hmat or smat",&
                                              hint=errmsg, calledby="eigen.F90")

eigen/eigen_hssetup.F90

@@ -61,7 +61,7 @@ CONTAINS
     DO i=1,nspins
        DO j=1,nspins
           CALL smat(i,j)%init(l_real,lapw%nv(i)+atoms%nlotot,lapw%nv(j)+atoms%nlotot,mpi%sub_comm,.false.)
-          CALL hmat(i,j)%init(l_real,lapw%nv(i)+atoms%nlotot,lapw%nv(j)+atoms%nlotot,mpi%sub_comm,.false.)
+          CALL hmat(i,j)%init(smat(i,j))
        ENDDO
     ENDDO
 
@@ -94,7 +94,7 @@ CONTAINS
     ! In collinear case only a copy is done
     ! In the parallel case also a redistribution happens
     CALL eigen_redist_matrix(mpi,lapw,atoms,smat,smat_final)
-    CALL eigen_redist_matrix(mpi,lapw,atoms,hmat,hmat_final)
+    CALL eigen_redist_matrix(mpi,lapw,atoms,hmat,hmat_final,smat_final)
     
   END SUBROUTINE eigen_hssetup
 END MODULE m_eigen_hssetup

eigen/eigen_redist_matrix.f90

@@ -14,7 +14,7 @@ CONTAINS
   !! In the non-collinear case, the 2x2 array of matrices is combined into the final matrix. Again a redistribution will happen in the parallel case
  
   
-  SUBROUTINE eigen_redist_matrix(mpi,lapw,atoms,mat,mat_final)
+  SUBROUTINE eigen_redist_matrix(mpi,lapw,atoms,mat,mat_final,mat_final_templ)
     USE m_types
     USE m_types_mpimat
     IMPLICIT NONE
@@ -23,33 +23,37 @@ CONTAINS
     TYPE(t_atoms),INTENT(IN)  :: atoms
     CLASS(t_mat),INTENT(INOUT):: mat(:,:)
     CLASS(t_mat),INTENT(INOUT):: mat_final
+    CLASS(t_mat),INTENT(IN),OPTIONAL :: mat_final_templ
 
     INTEGER:: m
 
     !determine final matrix size and allocate the final matrix
     m=lapw%nv(1)+atoms%nlotot
     IF (SIZE(mat)>1) m=m+lapw%nv(2)+atoms%nlotot
+    IF (.NOT.PRESENT(mat_final_templ)) THEN
        CALL mat_final%init(mat(1,1)%l_real,m,m,mpi%sub_comm,.TRUE.) !here the .true. creates a block-cyclic scalapack distribution
-    
+    ELSE
+       CALL mat_final%init(mat_final_templ)
+    ENDIF
     !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()
+       IF (.NOT.PRESENT(mat_final_templ)) CALL mat(1,1)%free()
        RETURN
     ENDIF
 
     CALL mat_final%copy(mat(1,1),1,1)
-    CALL mat(1,1)%free()
+    IF (.NOT.PRESENT(mat_final_templ)) CALL mat(1,1)%free()
   
     !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()
+    IF (.NOT.PRESENT(mat_final_templ)) CALL mat(2,2)%free()
 
     !Now collect off-diagonal parts
     CALL mat(1,2)%add_transpose(mat(2,1))
     CALL mat_final%copy(mat(1,2),1,lapw%nv(1)+atoms%nlotot+1)
-    CALL mat(1,2)%free()
-    CALL mat(2,1)%free()
+    IF (.NOT.PRESENT(mat_final_templ)) CALL mat(1,2)%free()
+    IF (.NOT.PRESENT(mat_final_templ)) CALL mat(2,1)%free()
     
   END SUBROUTINE eigen_redist_matrix
 END MODULE m_eigen_redist_matrix

types/types_mat.F90

@@ -24,7 +24,9 @@ MODULE m_types_mat
      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        :: init_details => t_mat_init
+     PROCEDURE        :: init_template => t_mat_init_template              !> initalize the matrix(overloaded for t_mpimat)
+     GENERIC          :: init => init_details,init_template
      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)
    END type t_mat
@@ -71,7 +73,28 @@ MODULE m_types_mat
 
      CALL mat%alloc(l_real,matsize1,matsize2)
    END SUBROUTINE t_mat_init
-
+   SUBROUTINE t_mat_init_template(mat,templ)
+     IMPLICIT NONE
+     CLASS(t_mat),INTENT(INOUT) :: mat
+     CLASS(t_mat),INTENT(IN)    :: templ
+     SELECT TYPE(templ)
+     TYPE is(t_mat)
+        mat%l_real=templ%l_real
+        mat%matsize1=templ%matsize1
+        mat%matsize2=templ%matsize2
+        IF (mat%l_real) THEN
+           ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2))
+           ALLOCATE(mat%data_c(1,1))
+           mat%data_r=0.0
+        ELSE
+           ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2))
+           ALLOCATE(mat%data_r(1,1))
+           mat%data_c=0.0
+        END IF
+     CLASS default
+        CALL judft_error("Mixed initialization in t_mat not possible(BUG)")
+     END SELECT
+   END SUBROUTINE t_mat_init_template
      
   SUBROUTINE t_mat_alloc(mat,l_real,matsize1,matsize2,init)
     CLASS(t_mat) :: mat

types/types_mpimat.F90

@@ -29,9 +29,9 @@ MODULE m_types_mpimat
      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   :: init_details => mpimat_init
+     PROCEDURE,PASS   :: init_template =>mpimat_init_template     !<overwriten from t_mat, also calls alloc in t_mat
      PROCEDURE,PASS   :: add_transpose => mpimat_add_transpose !<overwriten from t_mat
-
      PROCEDURE,PASS   :: generate_full_matrix    ! construct full matrix if only upper triangle of hermitian matrix is given
      PROCEDURE,PASS   :: print_matrix
      PROCEDURE,PASS   :: from_non_dist
@@ -286,6 +286,34 @@ CONTAINS
 #endif    
   END SUBROUTINE mpimat_init
 
+  SUBROUTINE mpimat_init_template(mat,templ)
+    IMPLICIT NONE
+    CLASS(t_mpimat),INTENT(INOUT)  :: mat
+    CLASS(t_mat),INTENT(IN)        :: templ
+
+    SELECT TYPE(templ)
+    TYPE IS (t_mpimat)
+       mat%l_real=templ%l_real
+       mat%matsize1=templ%matsize1
+       mat%matsize2=templ%matsize2
+       mat%global_size1=templ%global_size1
+       mat%global_size2=templ%global_size2
+       mat%mpi_com=templ%mpi_com
+       mat%blacs_desc=templ%blacs_desc
+       mat%blacs_ctext=templ%blacs_ctext
+       mat%npcol=templ%npcol
+       mat%nprow=templ%nprow
+       IF (mat%l_real) THEN
+          ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2))
+       ELSE
+          ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2))
+       END IF
+       CLASS default
+          CALL judft_error("Mixed initialization in t_mpimat not possible(BUG)")
+    END SELECT
+  END SUBROUTINE mpimat_init_template
+
+    
   SUBROUTINE priv_create_blacsgrid(mpi_subcom,l_2d,m1,m2,nbc,nbr,ictextblacs,sc_desc,local_size1,local_size2,nprow,npcol)
     IMPLICIT NONE
     INTEGER,INTENT(IN) :: mpi_subcom