Provide interfaces for PARPACK routines

vnquanvuong · Aug 20, 2024 · 97d4209 · 97d4209
1 parent 9401bea
commit 97d4209
Show file tree

Hide file tree

Showing 3 changed files with 166 additions and 94 deletions.
diff --git a/src/dftbp/extlibs/arpack.F90 b/src/dftbp/extlibs/arpack.F90
@@ -9,137 +9,217 @@
 
 #:assert not (WITH_ARPACK and INSTANCE_SAFE_BUILD)
 
+#:set ARPACK_PROC_PREFIX = "" if WITH_ARPACK else "module"
+#:set PARPACK_PROC_PREFIX = "" if WITH_PARPACK else "module"
+#:set PREFIXES_AND_KINDS =  [("s", "rsp"), ("d", "rdp")]
+
+
 !> Interfaces for the ARPACK routines needed in DFTB+ (currently for the linear response excited
 !> state calculations).
 module dftbp_extlibs_arpack
   use dftbp_common_accuracy, only : rsp, rdp
-#:if not WITH_ARPACK
-  use dftbp_io_message
-#:endif
   implicit none
-  private
 
+  private
   public :: withArpack, saupd, seupd
+  public :: withParpack, psaupd, pseupd
 
-#:if WITH_ARPACK
 
   !> Whether code was built with ARPACK support
-  logical, parameter :: withArpack = .true.
+  logical, parameter :: withArpack = ${FORTRAN_LOGICAL(WITH_ARPACK)}$
 
-#:else
+  !> Whether code was built with PARPACK support
+  logical, parameter :: withParpack = ${FORTRAN_LOGICAL(WITH_PARPACK)}$
 
-  !> Whether code was built with ARPACK support
-  logical, parameter :: withArpack = .false.
 
-  !> Dummy routines, as ARPACK library is not compiled in
+  ! Function overloading to be used within DFTB+
   interface saupd
-  #:for PREC in [("s"),("d")]
-    module procedure ${PREC}$saupd
+  #:for PREFIX, _ in PREFIXES_AND_KINDS
+    ${ARPACK_PROC_PREFIX}$ procedure ${PREFIX}$saupd
   #:endfor
-  end interface saupd
+  end interface
 
-  !> Dummy routines, as ARPACK library is not compiled in
+  ! Function overloading to be used within DFTB+
   interface seupd
-  #:for PREC in [("s"),("d")]
-    module procedure ${PREC}$seupd
+  #:for PREFIX, _ in PREFIXES_AND_KINDS
+    ${ARPACK_PROC_PREFIX}$ procedure ${PREFIX}$seupd
   #:endfor
   end interface
 
-contains
-
-  !> Generates error message, if a stub was called
-  subroutine stubError(routineName)
-    character(*), intent(in) :: routineName
-
-    call error("Internal error: " // trim(routineName) // "() called in a build without ARPACK&
-        & support")
+  ! Function overloading to be used within DFTB+
+  interface psaupd
+  #:for PREFIX, _ in PREFIXES_AND_KINDS
+    ${PARPACK_PROC_PREFIX}$ procedure p${PREFIX}$saupd
+  #:endfor
+  end interface
 
-  end subroutine stubError
+  ! Function overloading to be used within DFTB+
+  interface pseupd
+  #:for PREFIX, _ in PREFIXES_AND_KINDS
+    ${PARPACK_PROC_PREFIX}$ procedure p${PREFIX}$seupd
+  #:endfor
+  end interface
 
-#:endif
+  ! Interface definition of the routines
+  interface
+  #:for PREFIX, KIND in PREFIXES_AND_KINDS
 
-#:if WITH_ARPACK
-  !> Wrapper around ARPACK routines ssaupd/dsaupd.
-  interface saupd
-#:endif
-  #:for PREC, LABEL, VTYPE in [("s","single","rsp"),("d","double","rdp")]
-  #:if not WITH_ARPACK
-    !> Dummy ARPACK routine
-  #:endif
-    !> ${LABEL}$ precision Arnoldi solver call
-    subroutine ${PREC}$saupd(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,&
-        & workd, workl, lworkl, info)
+    !> Arnoldi solver call
+    ${ARPACK_PROC_PREFIX}$ subroutine ${PREFIX}$saupd(ido, bmat, n, which, nev, tol, resid, ncv,&
+        & v, ldv, iparam, ipntr, workd, workl, lworkl, info)
     #:if WITH_ARPACK
-      import :: ${VTYPE}$
+      import :: ${KIND}$
     #:endif
       integer, intent(inout) :: ido
       character, intent(in) :: bmat
       integer, intent(in) :: n
       character(2), intent(in) :: which
       integer, intent(in) :: nev
-      real(${VTYPE}$), intent(in) :: tol
-      real(${VTYPE}$), intent(inout) :: resid(n)
+      real(${KIND}$), intent(in) :: tol
+      real(${KIND}$), intent(inout) :: resid(n)
       integer, intent(in) :: ncv
       integer, intent(in) :: ldv
-      real(${VTYPE}$), intent(out) :: v(ldv, ncv)
+      real(${KIND}$), intent(out) :: v(ldv, ncv)
       integer, intent(inout) :: iparam(11)
       integer, intent(out) :: ipntr(11)
-      real(${VTYPE}$), intent(inout) :: workd(3 * n)
+      real(${KIND}$), intent(inout) :: workd(3 * n)
       integer, intent(in) :: lworkl
-      real(${VTYPE}$), intent(inout) :: workl(lworkl)
+      real(${KIND}$), intent(inout) :: workl(lworkl)
       integer, intent(inout) :: info
-     #:if not WITH_ARPACK
-      call stubError("${PREC}$saupd")
-     #:endif
-    end subroutine ${PREC}$saupd
-  #:endfor
-#:if WITH_ARPACK
-  end interface saupd
-#:endif
+    end subroutine ${PREFIX}$saupd
 
-#:if WITH_ARPACK
-  !> Wrapper around ARPACK routines sseupd/dseupd.
-  interface seupd
-#:endif
-  #:for PREC, LABEL, VTYPE in [("s","single","rsp"),("d","double","rdp")]
-  #:if not WITH_ARPACK
-    !> Dummy ARPACK routine
-  #:endif
-    !> ${LABEL}$ precision return from the results of the solver
-    subroutine ${PREC}$seupd(rvec, howmny, sel, d, z, ldz, sigma, bmat, n, which, nev, tol, resid,&
-        & ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info)
+    ${ARPACK_PROC_PREFIX}$ subroutine ${PREFIX}$seupd(rvec, howmny, sel, d, z, ldz, sigma, bmat,&
+        & n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info)
     #:if WITH_ARPACK
-      import :: ${VTYPE}$
+      import :: ${KIND}$
     #:endif
       logical, intent(in) :: rvec
       character, intent(in) :: howmny
       integer, intent(in) :: ncv
       logical, intent(in) :: sel(ncv)
       integer, intent(in) :: nev
-      real(${VTYPE}$), intent(out) :: d(nev)
+      real(${KIND}$), intent(out) :: d(nev)
       integer, intent(in) :: ldz
-      real(${VTYPE}$), intent(out) :: z(ldz, nev)
-      real(${VTYPE}$), intent(in) :: sigma
+      real(${KIND}$), intent(out) :: z(ldz, nev)
+      real(${KIND}$), intent(in) :: sigma
       character, intent(in) :: bmat
       integer, intent(in) :: n
       character(2), intent(in) :: which
-      real(${VTYPE}$), intent(in) :: tol
-      real(${VTYPE}$), intent(in) :: resid(n)
+      real(${KIND}$), intent(in) :: tol
+      real(${KIND}$), intent(in) :: resid(n)
       integer, intent(in) :: ldv
-      real(${VTYPE}$), intent(inout) :: v(ldv, ncv)
+      real(${KIND}$), intent(inout) :: v(ldv, ncv)
       integer, intent(in) :: iparam(7)
       integer, intent(inout) :: ipntr(11)
-      real(${VTYPE}$), intent(inout) :: workd(2 * n)
+      real(${KIND}$), intent(inout) :: workd(2 * n)
+      integer, intent(in) :: lworkl
+      real(${KIND}$), intent(inout) :: workl(lworkl)
+      integer, intent(inout) :: info
+    end subroutine ${PREFIX}$seupd
+
+    ${PARPACK_PROC_PREFIX}$ subroutine p${PREFIX}$saupd(comm, ido, bmat, n, which, nev, tol, resid,&
+        & ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info)
+    #:if WITH_PARPACK
+      import :: ${KIND}$
+    #:endif
+      integer, intent(in) :: comm
+      integer, intent(inout) :: ido
+      character, intent(in) :: bmat
+      integer, intent(in) :: n
+      character(2), intent(in) :: which
+      integer, intent(in) :: nev
+      real(${KIND}$), intent(in) :: tol
+      real(${KIND}$), intent(inout) :: resid(n)
+      integer, intent(in) :: ncv
+      integer, intent(in) :: ldv
+      real(${KIND}$), intent(out) :: v(ldv, ncv)
+      integer, intent(inout) :: iparam(11)
+      integer, intent(out) :: ipntr(11)
+      real(${KIND}$), intent(inout) :: workd(3 * n)
       integer, intent(in) :: lworkl
-      real(${VTYPE}$), intent(inout) :: workl(lworkl)
+      real(${KIND}$), intent(inout) :: workl(lworkl)
       integer, intent(inout) :: info
-     #:if not WITH_ARPACK
-      call stubError("${PREC}$seupd")
-     #:endif
-    end subroutine ${PREC}$seupd
+    end subroutine p${PREFIX}$saupd
+
+    ${PARPACK_PROC_PREFIX}$ subroutine p${PREFIX}$seupd(comm, rvec, howmny, sel, d, z, ldz, sigma,&
+        & bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info)
+    #:if WITH_PARPACK
+      import :: ${KIND}$
+    #:endif
+      integer, intent(in) :: comm
+      logical, intent(in) :: rvec
+      character, intent(in) :: howmny
+      integer, intent(in) :: ncv
+      logical, intent(inout) :: sel(ncv)
+      integer, intent(in) :: nev
+      real(${KIND}$), intent(out) :: d(nev)
+      integer, intent(in) :: ldz
+      real(${KIND}$), intent(out) :: z(ldz, nev)
+      real(${KIND}$), intent(in) :: sigma
+      character, intent(in) :: bmat
+      integer, intent(in) :: n
+      character(2), intent(in) :: which
+      real(${KIND}$), intent(in) :: tol
+      real(${KIND}$), intent(in) :: resid(n)
+      integer, intent(in) :: ldv
+      real(${KIND}$), intent(inout) :: v(ldv, ncv)
+      integer, intent(in) :: iparam(7)
+      integer, intent(inout) :: ipntr(11)
+      real(${KIND}$), intent(inout) :: workd(2 * n)
+      integer, intent(in) :: lworkl
+      real(${KIND}$), intent(inout) :: workl(lworkl)
+      integer, intent(inout) :: info
+    end subroutine p${PREFIX}$seupd
+
   #:endfor
-#:if WITH_ARPACK
-  end interface seupd
-#:endif
+  end interface
 
 end module dftbp_extlibs_arpack
+
+
+#:if (not WITH_ARPACK) or (not WITH_PARPACK)
+
+!> Defines stubs for ARPACK/PARPACK routines, in case the libraries are not present.
+submodule (dftbp_extlibs_arpack) dftbp_extlibs_arpack_stubs
+  use dftbp_io_message
+  implicit none
+
+contains
+
+#:for PREFIX, _ in PREFIXES_AND_KINDS
+
+#:if not WITH_ARPACK
+  module procedure ${PREFIX}$saupd
+    call stubError_("${PREFIX}$saupd", "ARPACK")
+  end procedure
+
+  module procedure ${PREFIX}$seupd
+    call stubError_("${PREFIX}$seupd", "ARPACK")
+  end procedure
+#:endif
+
+#:if not WITH_PARPACK
+  module procedure p${PREFIX}$saupd
+    call stubError_("p${PREFIX}$saupd", "PARPACK")
+  end procedure
+
+  module procedure p${PREFIX}$seupd
+    call stubError_("p${PREFIX}$seupd", "PARPACK")
+  end procedure
+#:endif
+
+#:endfor
+
+
+  !! Generates error message, if a stub was called
+  subroutine stubError_(routineName, libraryName)
+    character(*), intent(in) :: routineName, libraryName
+
+    call error("Internal error: " // routineName // "() called in a build without " // libraryName&
+        & // " suppport")
+
+  end subroutine stubError_
+
+end submodule dftbp_extlibs_arpack_stubs
+
+#:endif
diff --git a/src/dftbp/include/common.fypp b/src/dftbp/include/common.fypp
@@ -24,6 +24,7 @@
 #:set WITH_MBD = defined('WITH_MBD')
 #:set WITH_SOCKETS = defined('WITH_SOCKETS')
 #:set WITH_ARPACK = defined('WITH_ARPACK')
+#:set WITH_PARPACK = WITH_ARPACK and WITH_MPI
 #:set WITH_PLUMED = defined('WITH_PLUMED')
 #:set WITH_TRANSPORT = defined('WITH_TRANSPORT')
 #:set WITH_POISSON = defined('WITH_POISSON')

diff --git a/src/dftbp/timedep/linrespgrad.F90 b/src/dftbp/timedep/linrespgrad.F90
@@ -23,7 +23,7 @@ module dftbp_timedep_linrespgrad
   use dftbp_dftb_shortgammafuncs, only : expGammaPrime
   use dftbp_dftb_sk, only : rotateH0
   use dftbp_dftb_slakocont, only : TSlakoCont, getMIntegrals, getSKIntegrals
-  use dftbp_extlibs_arpack, only : withArpack, saupd, seupd
+  use dftbp_extlibs_arpack, only : psaupd, pseupd, saupd, seupd, withArpack
   use dftbp_io_message, only : error
   use dftbp_io_taggedoutput, only : TTaggedWriter, tagLabels
   use dftbp_math_blasroutines, only : gemm, hemv, symm, herk
@@ -970,8 +970,7 @@ subroutine buildAndDiagExcMatrixArpack(iGlobal, fGlobal, env, orb, lr, rpa, tran
     character(lc) :: tmpStr
     type(TFileDescr) :: fdArnoldiTest
 
-  #:if WITH_SCALAPACK and WITH_ARPACK 
-    external pdsaupd, pdseupd
+  #:if WITH_PARPACK
     comm = env%mpi%globalComm%id 
   #:endif
 
@@ -1010,20 +1009,16 @@ subroutine buildAndDiagExcMatrixArpack(iGlobal, fGlobal, env, orb, lr, rpa, tran
     iparam(3) = maxArIter
     ! solve A*x = lambda*x, with A symmetric
     iparam(7) = 1
-    
+
     do
 
       ! call the reverse communication interface from arpack
-    #:if WITH_SCALAPACK and WITH_ARPACK
-
-      call pdsaupd(comm, ido, "I", nLoc, "SM", nexc, arTol, resid, ncv, vv, nLoc, iparam,&
+    #:if WITH_PARPACK
+      call psaupd(comm, ido, "I", nLoc, "SM", nexc, arTol, resid, ncv, vv, nLoc, iparam,&
           & ipntr, workd, workl, lworkl, info)
-
     #:else
-
       call saupd(ido, "I", rpa%nxov_rd, "SM", nexc, arTol, resid, ncv, vv, rpa%nxov_rd, iparam,&
           & ipntr, workd, workl, lworkl, info)
-
     #:endif
 
       if (ido == 99) then
@@ -1061,9 +1056,9 @@ subroutine buildAndDiagExcMatrixArpack(iGlobal, fGlobal, env, orb, lr, rpa, tran
       ! to DSAUPD.  These arguments MUST NOT BE MODIFIED between the the last call to DSAUPD and the
       ! call to DSEUPD.s
       ! Note: At this point xpy holds the hermitian eigenvectors F
-    #:if WITH_SCALAPACK and WITH_ARPACK
+    #:if WITH_PARPACK
 
-      call pdseupd (comm, rvec, "All", selection, eval, vv, nLoc, sigma, "I", nLoc,& 
+      call pseupd (comm, rvec, "All", selection, eval, vv, nLoc, sigma, "I", nLoc,&
           & "SM", nexc, arTol, resid, ncv, vv, nLoc, iparam, ipntr, workd, workl, lworkl, info)
 
       xpy(:,:) = 0.0_dp
@@ -1212,10 +1207,6 @@ subroutine buildAndDiagExcMatrixStratmann(iGlobal, fGlobal, env, orb, lr, rpa, t
 
     logical :: didConverge
 
-  #:if WITH_SCALAPACK
-    external pdsaupd, pdseupd
-  #:endif
-
     ! Local chunk of RPA vectors have this size under MPI
     nLoc = fGlobal - iGlobal + 1