GPU solver for modes

Co-authored-by: Tammo van der Heide <[email protected]>

CHANGELOG.rst

@@ -33,6 +33,8 @@ Added
 - More control over output of data and band structures during MD
   calculations
 
+- MAGMA GPU accelerated solver for the modes code
+
 Fixed
 -----
 

app/modes/initmodes.F90

@@ -14,6 +14,9 @@ module modes_initmodes
   use dftbp_common_file, only : TFileDescr, openFile, closeFile
   use dftbp_common_filesystem, only : findFile, joinPathsPrettyErr, getParamSearchPaths
   use dftbp_common_globalenv, only : stdOut
+#:if WITH_MAGMA
+  use dftbp_common_gpuenv, only : TGpuEnv, TGpuEnv_init
+#:endif
   use dftbp_common_release, only : releaseYear
   use dftbp_common_unitconversion, only : massUnits
   use dftbp_extlibs_xmlf90, only : fnode, fNodeList, string, char, getLength, getItem1,&
@@ -115,13 +118,21 @@ module modes_initmodes
   !> Eigensolver choice
   integer :: iSolver
 
+#:if WITH_MAGMA
+  !> Global GPU settings
+  type(TGpuEnv), public :: gpu
+#:endif
+
   !> Namespace for possible eigensolver methods
   type :: TSolverTypesEnum
 
     ! lapack/scalapack solvers
     integer :: qr = 1
-    integer :: divideandconquer = 2
-    integer :: relativelyrobust = 3
+    integer :: divideAndConquer = 2
+    integer :: relativelyRobust = 3
+
+    ! GPU accelerated solver using MAGMA
+    integer :: magmaEvd = 4
 
   end type TSolverTypesEnum
 
@@ -204,9 +215,16 @@ subroutine initProgramVariables()
     case ("qr")
       iSolver = solverTypes%qr
     case ("divideandconquer")
-      iSolver = solverTypes%divideandconquer
+      iSolver = solverTypes%divideAndConquer
     case ("relativelyrobust")
-      iSolver = solverTypes%relativelyrobust
+      iSolver = solverTypes%relativelyRobust
+    case ("magma")
+    #:if WITH_MAGMA
+      call TGpuEnv_init(gpu)
+    #:else
+      call error("Magma-solver selected, but program was compiled without MAGMA")
+    #:endif
+      iSolver = solverTypes%magmaEvd
     case default
       call detailedError(root, "Unknown eigensolver "//char(buffer2))
     end select

app/modes/modes.F90

@@ -17,10 +17,16 @@ program modes
   use dftbp_io_message, only : error
   use dftbp_io_taggedoutput, only : TTaggedWriter, TTaggedWriter_init
   use dftbp_math_eigensolver, only : heev, heevd, heevr
+#:if WITH_MAGMA
+  use dftbp_math_eigensolver, only : gpu_evd
+#:endif
   use modes_initmodes, only : dynMatrix, bornMatrix, bornDerivsMatrix, modesToPlot, geo,&
       & iMovedAtoms, nCycles, nDerivs, nModesToPlot, nMovedAtom, nSteps, tAnimateModes, tPlotModes,&
       & tEigenVectors, tRemoveRotate, tRemoveTranslate, atomicMasses, initProgramVariables,&
       & iSolver, solverTypes, setEigvecGauge
+#:if WITH_MAGMA
+  use modes_initmodes, only : gpu
+#:endif
   use modes_modeprojection, only : project
 #:if WITH_MPI
   use mpi, only : MPI_THREAD_FUNNELED
@@ -92,10 +98,16 @@ program modes
   select case(iSolver)
   case(solverTypes%qr)
     call heev(dynMatrix, eigenValues, "U", eigenSolverMode)
-  case(solverTypes%divideandconquer)
+  case(solverTypes%divideAndConquer)
     call heevd(dynMatrix, eigenValues, "U", eigenSolverMode)
-  case(solverTypes%relativelyrobust)
+  case(solverTypes%relativelyRobust)
     call heevr(dynMatrix, eigenValues, "U", eigenSolverMode)
+  case(solverTypes%magmaEvd)
+  #:if WITH_MAGMA
+    call gpu_evd(gpu%nGpu, dynMatrix, eigenValues, "U", eigenSolverMode)
+  #:else
+    call error("Magma-solver selected, but program was compiled without MAGMA")
+  #:endif
   case default
     call error("Unknown eigensolver choice")
   end select

doc/dftb+/manual/modes.tex

@@ -101,7 +101,9 @@ \subsubsection{DirectRead\cb}
   Slater-Koster files. See p.~\pref{sec:modes.Masses}.
 \item[\is{EigenSolver}] Choice of solver for the dynamical matrix,
   current choices are \kw{QR}, \kw{DivideAndConquer} or
-  \kw{RelativelyRobust}.
+  \kw{RelativelyRobust}. If the modes code is compiled with the
+  MAGMA~\cite{tdb10, tnld10, dghklty14} library
+  included, then the \kw{MAGMA} solver is also available.
 \end{description}
 
 

src/dftbp/elecsolvers/elecsolvertypes.F90

@@ -34,7 +34,7 @@ module dftbp_elecsolvers_elecsolvertypes
     integer :: gf = 11
     integer :: onlyTransport = 12
 
-    ! GPU accellerated solvers using MAGMA
+    ! GPU accelerated solvers using MAGMA
     integer :: magma_gvd = 13
 
   end type TElecSolverTypesEnum

src/dftbp/extlibs/magma.F90

@@ -11,7 +11,8 @@
 module dftbp_extlibs_magma
   use, intrinsic :: iso_c_binding, only : c_int
 #:if WITH_MAGMA
-  use magma, only : magmaf_ssygvd_m, magmaf_dsygvd_m, magmaf_chegvd_m, magmaf_zhegvd_m
+  use magma, only : magmaf_ssygvd_m, magmaf_dsygvd_m, magmaf_chegvd_m, magmaf_zhegvd_m,&
+      & magmaf_ssyevd_m, magmaf_dsyevd_m, magmaf_cheevd_m, magmaf_zheevd_m
 #:endif
   implicit none
 
@@ -20,6 +21,7 @@ module dftbp_extlibs_magma
 #:if WITH_MAGMA
   public :: getGpusAvailable, getGpusRequested, gpusInit
   public :: magmaf_ssygvd_m, magmaf_dsygvd_m, magmaf_chegvd_m, magmaf_zhegvd_m
+  public :: magmaf_ssyevd_m, magmaf_dsyevd_m, magmaf_cheevd_m, magmaf_zheevd_m
 #:endif
 
   !> Whether code was built with GPU support

src/dftbp/math/eigensolver.F90

@@ -19,14 +19,15 @@ module dftbp_math_eigensolver
   use dftbp_io_message, only : error, warning
 #:if WITH_MAGMA
   use dftbp_extlibs_magma,  only : magmaf_ssygvd_m, magmaf_dsygvd_m, magmaf_chegvd_m,&
-      & magmaf_zhegvd_m
+      & magmaf_zhegvd_m, magmaf_ssyevd_m, magmaf_dsyevd_m, magmaf_cheevd_m,&
+      & magmaf_zheevd_m
 #:endif
   implicit none
 
   private
   public :: heev, heevd, heevr, hegv, hegvd, gvr, geev
 #:if WITH_MAGMA
-  public :: gpu_gvd
+  public :: gpu_gvd, gpu_evd
 #:endif
 
 
@@ -96,13 +97,23 @@ module dftbp_math_eigensolver
 
 
 #:if WITH_MAGMA
-  !> Divide and conquer MAGMA GPU eigensolver
+
+  !> Divide and conquer MAGMA GPU generalised eigensolver
   interface gpu_gvd
     module procedure real_magma_ssygvd
     module procedure dble_magma_dsygvd
     module procedure cmplx_magma_chegvd
     module procedure dblecmplx_magma_zhegvd
-  end interface
+  end interface gpu_gvd
+
+  !> Divide and conquer MAGMA GPU eigensolver
+  interface gpu_evd
+    module procedure real_magma_ssyevd
+    module procedure dble_magma_dsyevd
+    module procedure cmplx_magma_cheevd
+    module procedure dblecmplx_magma_zheevd
+  end interface gpu_evd
+
 #:endif
 
   !> Simple eigensolver for a general matrix
@@ -1985,6 +1996,105 @@ end subroutine dblecmplx_zhegvr
 
   end subroutine ${DTYPE}$_magma_${NAME}$
 
+#:endfor
+
+  #:for DTYPE, VPREC, VTYPE, NAME in [('real', 's', 'real', 'ssyevd'),&
+  & ('dble', 'd', 'real', 'dsyevd'), ('cmplx', 's', 'complex', 'cheevd'),&
+  & ('dblecmplx', 'd', 'complex', 'zheevd')]
+  !> Eigensolution for symmetric/hermitian matrices on GPU(s)
+  subroutine ${DTYPE}$_magma_${NAME}$(ngpus, a, w, uplo, jobz)
+
+    !> Number of GPUs to use
+    integer, intent(in) :: ngpus
+
+    !> contains the matrix for the solver, returns eigenvectors if requested (matrix always
+    !! overwritten on return anyway)
+    ${VTYPE}$(r${VPREC}$p), intent(inout) :: a(:,:)
+
+    !> eigenvalues
+    real(r${VPREC}$p), intent(out) :: w(:)
+
+    !> upper or lower triangle of the matrix
+    character, intent(in) :: uplo
+
+    !> compute eigenvalues 'N' or eigenvalues and eigenvectors 'V'
+    character, intent(in) :: jobz
+
+    ${VTYPE}$(r${VPREC}$p), allocatable :: work(:)
+
+  #:if VTYPE == 'complex'
+    real(r${VPREC}$p), allocatable :: rwork(:)
+    integer :: lrwork
+  #:endif
+
+    integer, allocatable :: iwork(:)
+    integer :: lwork, liwork, n, info
+    character(len=100) :: error_string
+
+    @:ASSERT(uplo == 'u' .or. uplo == 'U' .or. uplo == 'l' .or. uplo == 'L')
+    @:ASSERT(jobz == 'n' .or. jobz == 'N' .or. jobz == 'v' .or. jobz == 'V')
+    @:ASSERT(all(shape(a)==size(w)))
+    n = size(a, dim=1)
+    @:ASSERT(n>0)
+
+    ! Workspace query
+    allocate(work(1))
+    allocate(iwork(1))
+  #:if VTYPE == 'complex'
+    allocate(rwork(1))
+  #:endif
+
+    call magmaf_${NAME}$_m(ngpus, jobz, uplo, n, a, n, w, work, -1,&
+      #:if VTYPE == 'complex'
+        & rwork, -1,&
+      #:endif
+        & iwork, -1, info)
+
+    if (info /= 0) then
+      call error("Failure in MAGMA_${NAME}$ to determine optimum workspace")
+    endif
+
+ #:if VTYPE == 'complex'
+    lwork = floor(real(work(1)))
+    lrwork = floor(rwork(1))
+    liwork = int(iwork(1))
+    deallocate(work) ;  allocate(work(lwork))
+    deallocate(rwork) ;  allocate(rwork(lrwork))
+    deallocate(iwork) ; allocate(iwork(liwork))
+  #:endif
+  #:if VTYPE == 'real'
+    lwork = floor(work(1))
+    liwork = int(iwork(1))
+    deallocate(work) ;  allocate(work(lwork))
+    deallocate(iwork) ; allocate(iwork(liwork))
+   #:endif
+
+    ! MAGMA diagonalization
+    call magmaf_${NAME}$_m(ngpus, jobz, uplo, n, a, n, w, work, lwork,&
+      #:if VTYPE == 'complex'
+        & rwork, lrwork,&
+      #:endif
+        & iwork, liwork, info)
+
+    ! test for errors
+    if (info /= 0) then
+      if (info < 0) then
+        write(error_string, "('Failure in diagonalisation routine magmaf_${NAME}$_m,&
+            & illegal argument at position ',i6)") info
+        call error(error_string)
+      else if (info <= n) then
+        write(error_string, "('Failure in diagonalisation routine magmaf_${NAME}$_m,&
+            & diagonal element ',i6,' did not converge to zero.')") info
+        call error(error_string)
+      else
+        write(error_string, "('Failure in diagonalisation routine magmaf_${NAME}$_m,&
+            & non-positive definite overlap! ',i6)") info - n
+        call error(error_string)
+      endif
+    endif
+
+  end subroutine ${DTYPE}$_magma_${NAME}$
+
 #:endfor
 
 #:endif

test/app/modes/C6H6_gpu/_vibrations.tag

@@ -0,0 +1,13 @@
+frequencies         :real:1:36
+ -0.105408819839487E-009 -0.613898933299534E-010 -0.506005818810573E-010
+  0.464266256879880E-010  0.675455232232590E-010  0.103727588699306E-009
+  0.170796258550142E-002  0.170796274712918E-002  0.283235854111502E-002
+  0.283235861655495E-002  0.297515692345838E-002  0.310763510055179E-002
+  0.376571149581518E-002  0.376571164383678E-002  0.410515844274293E-002
+  0.412386823863947E-002  0.412386837443244E-002  0.467555768174623E-002
+  0.501536479599792E-002  0.501536494343011E-002  0.523030635218502E-002
+  0.525155968206473E-002  0.531224661080881E-002  0.531224671591227E-002
+  0.608112043290847E-002  0.678193393596265E-002  0.723326049312746E-002
+  0.723326057601041E-002  0.827283834993114E-002  0.827283841047394E-002
+  0.139562576095333E-001  0.139739592302981E-001  0.139739593600197E-001
+  0.140256517237074E-001  0.140256518626061E-001  0.140563193140046E-001

test/app/modes/C6H6_gpu/geom.gen

@@ -0,0 +1,14 @@
+ 12  C
+ C H
+    1  1   -0.2906476665E-12    0.4177079156E+00    0.1341743646E+01
+    2  1   -0.5269647781E-12    0.1364195076E+01    0.3129614748E+00
+    3  1   -0.5720571923E-12    0.9464871600E+00   -0.1021111535E+01
+    4  1   -0.4020479748E-12   -0.4177079139E+00   -0.1326402376E+01
+    5  1    0.1721617016E-12   -0.1364195075E+01   -0.2976202048E+00
+    6  1    0.2335252272E-12   -0.9464871584E+00    0.1036452805E+01
+    7  2    0.7860404010E-12    0.7453789618E+00    0.2388257444E+01
+    8  2    0.6582302911E-12    0.2434338133E+01    0.5524469271E+00
+    9  2    0.7307045197E-13    0.1688959174E+01   -0.1828139883E+01
+   10  2    0.7006607460E-13   -0.7453789602E+00   -0.2372916174E+01
+   11  2   -0.5559460858E-12   -0.2434338131E+01   -0.5371056572E+00
+   12  2    0.3545695500E-12   -0.1688959172E+01    0.1843481153E+01