Bugfixes

diagonalization/CMakeLists.txt

@@ -7,7 +7,9 @@ diagonalization/magma.F90
 diagonalization/scalapack.F90
 diagonalization/chase_diag.F90
 diagonalization/symmetrize_matrix.f90
+diagonalization/cusolver_diag.F90
 diagonalization/elemental.F90)
+
 if (FLEUR_USE_ELPA_20180525)
  set(fleur_F90 ${fleur_F90}
   diagonalization/elpa_20180525.F90
@@ -17,8 +19,4 @@ else()
   diagonalization/elpa.F90
 )
 endif()
-if (FLEUR_USE_GPU)
-set(fleur_F90 ${fleur_F90}
-  diagonalization/cusolver_diag.F90
-)
-endif()
+

diagonalization/cusolver_diag.F90

@@ -39,6 +39,7 @@ MODULE m_cusolver_diag
 CONTAINS
   SUBROUTINE cusolver_diag(hmat,smat,ne,eig,zmat)
     !Simple driver to solve Generalized Eigenvalue Problem using CuSolverDN
+    USE m_types
     IMPLICIT NONE
     CLASS(t_mat),INTENT(INOUT) :: hmat,smat
     INTEGER,INTENT(INOUT)      :: ne
@@ -49,27 +50,14 @@ CONTAINS
     INTEGER,PARAMETER:: max_sweeps=15
     REAL             :: tol=1E-7
     
-    TYPE(t_gpumat)::smat_gpu,hmat_gpu
     
     ALLOCATE(t_mat::zmat)
     CALL zmat%alloc(hmat%l_real,hmat%matsize1,ne)
-    SELECT TYPE(hmat)
-    TYPE IS (t_mat)
-       CALL hmat_gpu%init_from(hmat)
-       CALL smat_gpu%init_from(smat)
-       IF (hmat%l_real) THEN
-          CALL cusolver_real(hmat_gpu%gpu_r,smat_gpu%gpu_r,smat%matsize1,ne,tol,max_sweeps,eig,z%data_r)
-       ELSE
-          CALL cusolver_complex(hmat_gpu%gpu_c,smat_gpu%gpu_c,smat%matsize1,ne,tol,max_sweeps,eig,z%data_c)
-       END IF
-    TYPE IS (t_gpumat)
-       IF (hmat%l_real) THEN
-          CALL cusolver_real(hmat%gpu_r,smat%gpu_r,smat%matsize1,ne,tol,max_sweeps,eig,z%data_r)
-       ELSE
-          CALL cusolver_complex(hmat%gpu_c,smat%gpu_c,smat%matsize1,ne,tol,max_sweeps,eig,z%data_c)
-       END IF
-       
-    END SELECT
+    IF (hmat%l_real) THEN
+       CALL cusolver_real(hmat%data_r,smat%data_r,smat%matsize1,ne,tol,max_sweeps,eig,zmat%data_r)
+    ELSE
+       CALL cusolver_complex(hmat%data_c,smat%data_c,smat%matsize1,ne,tol,max_sweeps,eig,zmat%data_c)
+    END IF
 #endif
        
   END SUBROUTINE cusolver_diag

diagonalization/eigen_diag.F90

@@ -29,8 +29,18 @@ MODULE m_eigen_diag
 #else
   INTEGER,PARAMETER:: diag_magma=-6
 #endif
-  INTEGER,PARAMETER:: diag_lapack=4
+#ifdef CPP_CHASE
   INTEGER,PARAMETER:: diag_chase=7
+#else
+  INTEGER,PARAMETER:: diag_chase=-7
+#endif
+#ifdef CPP_GPU
+  INTEGER,PARAMETER:: diag_cusolver=8
+#else
+  INTEGER,PARAMETER:: diag_cusolver=-8
+#endif
+  
+  INTEGER,PARAMETER:: diag_lapack=4
   INTEGER,PARAMETER:: diag_debugout=99
   PUBLIC eigen_diag,parallel_solver_available
 CONTAINS
@@ -47,7 +57,9 @@ CONTAINS
     USE m_elemental
     USE m_chase_diag
     USE m_types_mpimat
+    USE m_types_gpumat
     USE m_matrix_copy
+    USE m_cusolver_diag
     IMPLICIT NONE
 #ifdef CPP_MPI
     include 'mpif.h'
@@ -99,14 +111,12 @@ CONTAINS
        CALL scalapack(hmat,smat,ne,eig,ev)
     CASE (diag_magma)
        !CALL magma_diag(hmat,smat,ne,eig,ev)
+    CASE (diag_cusolver)
+       CALL cusolver_diag(hmat,smat,ne,eig,ev)
     CASE (diag_lapack)
        CALL lapack_diag(hmat,smat,ne,eig,ev)
     CASE (diag_chase)
-#ifdef CPP_CHASE
        CALL chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,ev)
-#else
-       CALL juDFT_error('ChASE eigensolver selected but not available', calledby = 'eigen_diag')
-#endif
     CASE (diag_debugout)
        CALL priv_debug_out(smat,hmat)
     END SELECT
@@ -215,11 +225,7 @@ CONTAINS
        diag_solver=diag_scalapack
 #endif
     ELSE
-#ifdef CPP_MAGMA
-       diag_solver=diag_magma
-#else
        diag_solver=diag_lapack
-#endif
     ENDIF
 
     !check if a special solver was requested
@@ -229,49 +235,17 @@ CONTAINS
     IF (trim(juDFT_string_for_argument("-diag"))=="lapack")    diag_solver=diag_lapack
     IF (trim(juDFT_string_for_argument("-diag"))=="magma")     diag_solver=diag_magma
     IF (trim(juDFT_string_for_argument("-diag"))=="chase")     diag_solver=diag_chase
+    IF (trim(juDFT_string_for_argument("-diag"))=="cusolver")  diag_solver=diag_cusolver
     IF (trim(juDFT_string_for_argument("-diag"))=="debugout")  diag_solver=diag_debugout
     
     !Check if solver is possible
-    if (diag_solver<0) call priv_solver_error(diag_solver,parallel)
-    IF (ANY((/diag_lapack,diag_magma/)==diag_solver).AND.parallel) CALL priv_solver_error(diag_solver,parallel)
-    if (any((/diag_elpa,diag_elemental,diag_scalapack/)==diag_solver).and..not.parallel) call priv_solver_error(diag_solver,parallel)
+    IF (diag_solver<0)  CALL juDFT_error("You selected a solver for the eigenvalue problem that is not available",hint="You most probably did not provide the appropriate libraries for compilation/linking")
+    IF (ANY((/diag_lapack,diag_magma,diag_cusolver/)==diag_solver).AND.parallel) CALL judft_error("You selected an eigensolver that does not support distributed memory parallism",hint="Try scalapack,elpa or another supported solver for parallel matrices")
+    IF (ANY((/diag_elpa,diag_elemental,diag_scalapack/)==diag_solver).AND..NOT.parallel) CALL judft_error("You selected an eigensolver for matrices that are memory distributed",hint="Try lapack, cusolver or another supported solver for non-distributed matrices")
 
   END FUNCTION priv_select_solver
 
 
-  SUBROUTINE priv_solver_error(diag_solver,parallel)
-    IMPLICIT NONE
-    INTEGER,INTENT(IN):: diag_solver
-    LOGICAL,INTENT(IN)::parallel
-    SELECT CASE(diag_solver)
-    CASE (diag_elpa)
-       IF (parallel) THEN
-          CALL juDFT_error("You did not compile with the ELPA solver and hence can not use it")
-       ELSE
-          CALL juDFT_error("The ELPA solver can not be used in serial")
-       ENDIF
-    CASE (diag_elemental)
-       IF (parallel) THEN
-          CALL juDFT_error("You did not compile with the ELEMENTAL solver and hence can not use it")
-       ELSE
-          CALL juDFT_error("The ELEMENTAL solver can not be used in serial")
-       ENDIF
-    CASE (diag_scalapack)
-       IF (parallel) THEN
-          CALL juDFT_error("You did not compile with the SCALAPACK solver and hence can not use it")
-       ELSE
-          CALL juDFT_error("The SCALAPACK solver can not be used in serial")
-       ENDIF
-    CASE (diag_lapack)
-       IF (parallel) THEN
-          CALL juDFT_error("The LAPACK solver can not be used in parallel")
-       ENDIF
-    CASE (diag_magma)
-       CALL juDFT_error("You have not compiled with MAGMA support")
-    CASE DEFAULT
-       CALL judft_error("You have selected an unkown eigensolver")
-    END SELECT
-  END SUBROUTINE priv_solver_error
-
+ 
 
 END MODULE m_eigen_diag

eigen/eigen_hssetup.F90

@@ -19,6 +19,7 @@ CONTAINS
        stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat_final,hmat_final)
     USE m_types
     USE m_types_mpimat
+    USE m_types_gpumat
     USE m_hs_int
     USE m_hsvac
     USE m_od_hsvac
@@ -92,7 +93,8 @@ CONTAINS
     ! Collect the four noco parts into a single matrix
     ! In collinear case only a copy is done
     ! In the parallel case also a redistribution happens
-    ALLOCATE(smat_final,hmat_final,source=smat(1,1))
+    ALLOCATE(smat_final,source=smat(1,1))
+    ALLOCATE(hmat_final,source=smat(1,1))
     CALL eigen_redist_matrix(mpi,lapw,atoms,smat,smat_final)
     CALL eigen_redist_matrix(mpi,lapw,atoms,hmat,hmat_final,smat_final)
     

types/CMakeLists.txt

@@ -24,6 +24,7 @@ types/types_regionCharges.f90
 types/types_dos.f90
 types/types_denCoeffsOffdiag.f90
 types/types_force.f90
+types/types_gpumat.F90
 )
 
 set(inpgen_F90 ${inpgen_F90}
@@ -50,8 +51,3 @@ types/types_dos.f90
 types/types_denCoeffsOffdiag.f90
 types/types_force.f90
 )
-if (FLEUR_USE_GPU)
- set(fleur_F90 ${fleur_F90}
- types/types_gpumat.F90
- )
-endif()	

types/types_gpumat.F90

@@ -7,6 +7,7 @@ MODULE m_types_gpumat
   !!
   !! Several select type constructs are only needed to make the PGI compiler happy...
   TYPE,EXTENDS(t_mat) :: t_gpumat
+#ifdef CPP_GPU     
    CONTAINS
      PROCEDURE        :: alloc =>t_gpumat_alloc
      PROCEDURE        :: multiply=>t_gpumat_multiply            !> do a matrix-matrix multiply
@@ -19,10 +20,11 @@ MODULE m_types_gpumat
      PROCEDURE        :: move => t_gpumat_move                  !> move data into another t_mat 
      PROCEDURE        :: add_transpose => t_gpumat_add_transpose!> add the tranpose/Hermitian conjg. without the diagonal 
      PROCEDURE        :: init_from                              !> Initialize from an existing matrix
+#endif
   END TYPE t_gpumat
 
 CONTAINS
-
+#ifdef CPP_GPU
   SUBROUTINE t_gpumat_alloc(mat,l_real,matsize1,matsize2,init)
     CLASS(t_gpumat) :: mat
     LOGICAL,INTENT(IN),OPTIONAL:: l_real
@@ -294,4 +296,5 @@ CONTAINS
        ENDIF
     END SELECT
   END SUBROUTINE t_gpumat_clear
+#endif  
 END MODULE m_types_gpumat