subroutine ylm4 ported onto GPU

cmake/Files_and_Targets.txt

@@ -56,7 +56,7 @@ kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F glo
  )
 
 set(inpgen_F90  ${inpgen_F90} global/constants.f90 io/xsf_io.f90
-eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.f90
+eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.F90
 global/enpara.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
 inpgen/closure.f90 inpgen/inpgen_arguments.F90
 juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 kpoints/julia.f90 global/utility.F90 

math/CMakeLists.txt

@@ -31,7 +31,7 @@ math/SphBessel.f90
 math/DoubleFactorial.f90
 math/ExpSave.f90
 math/intgr.F90
-math/ylm4.f90
+math/ylm4.F90
 )
 if (FLEUR_USE_FFTMKL)
    set(fleur_F90 ${fleur_F90} math/mkl_dfti.f90)

math/ylm4.f90 → math/ylm4.F90

       MODULE m_ylm
         IMPLICIT NONE
-        private
+        PRIVATE
 !---> check whether  or not normalizations are needed
         REAL, ALLOCATABLE, SAVE  :: ynorm(:)
         INTEGER,           SAVE  :: lmaxd = -1  ! initial value
-        public ylm4,ylmnorm_init
+        PUBLIC ylm4,ylmnorm_init
+#ifdef _CUDA
+        PUBLIC ylm4_dev
+#endif
       CONTAINS 
+
+#ifdef _CUDA
+      ATTRIBUTES(device) SUBROUTINE ylm4_dev(lmax,v,ylm)
+!************************************************************
+!     generate the spherical harmonics for the vector v
+!     using a stable upward recursion in l.  (see notes
+!     by m. weinert.)
+!          m.weinert   january 1982
+!     modified by R. Podloucky (added in ynorm); July 1989
+!     cleaned up    mw 1995
+!
+!     modified to make use of f90 constructs. note that
+!     the normalization is an internal subroutine and hence
+!     can only be called from here. also, no need to dimension
+!     arrays for ynorm, done dynamically.          mw 1999
+!************************************************************
+!      USE m_juDFT
+      INTEGER, VALUE, INTENT (IN) :: lmax
+      REAL,    DEVICE, INTENT (IN) :: v(3)
+      COMPLEX, DEVICE, INTENT (OUT):: ylm( (lmax+1)**2 )
+
+      REAL, PARAMETER   :: small = 1.0e-12
+
+      INTEGER  :: l,lm0,m
+      REAL     :: fac,x,y,z,xy,r,rxy,cth,sth,cph,sph,cph2
+      REAL,allocatable,device     :: c(:),s(:)
+      REAL,allocatable,device     :: p(:,:)
+      COMPLEX  :: ylms
+      REAL,allocatable,device     :: ynorm_dev(:)
+      REAL    a,cd,fpi
+
+      allocate(c(0:lmax),s(0:lmax),ynorm_dev((lmax+1)**2))
+      allocate(p(0:lmax,0:lmax))
+
+!--->    calculate norm
+
+      fpi = 4.0 * 3.1415926535897932  
+      DO l=0,lmax
+         lm0 = l*(l+1) + 1
+         cd=1.0
+         a = sqrt( (2*l+1)/fpi )
+         ynorm_dev(lm0) = a
+         DO m=1,l
+            cd=cd/( (l+1-m)*(l+m) )
+            ynorm_dev(lm0+m)  = a*sqrt(cd)
+            ynorm_dev(lm0-m) = ( (-1.0)**m )*ynorm_dev(lm0+m)
+         ENDDO
+      ENDDO
+
+!--->    calculate sin and cos of theta and phi
+      x = v(1)
+      y = v(2)
+      z = v(3)
+      xy  = x*x + y*y
+      r   = sqrt(xy + z*z)
+      rxy = sqrt(xy)
+
+      IF (r.GT.small) THEN
+         cth = z/r
+         sth = rxy/r
+      ELSE
+         sth = 0.0
+         cth = 1.0
+      ENDIF
+      IF(rxy.GT.small) THEN
+         cph = x/rxy
+         sph = y/rxy
+      ELSE
+         cph = 1.0
+         sph = 0.0
+      ENDIF
+
+!---> generate associated legendre functions for m.ge.0
+      fac = 1.0
+!---> loop over m values
+      DO m = 0, lmax - 1
+         fac = -(m+m-1)*fac
+         p(m,m) = fac
+         p(m+1,m) = (m+m+1)*cth*fac
+!--->    recurse upward in l
+         DO l = m+2, lmax
+            p(l,m)=((l+l-1)*cth*p(l-1,m)-(l+m-1)*p(l-2,m))/(l-m)
+         ENDDO
+         fac = fac*sth
+      ENDDO
+      p(lmax,lmax) = -(lmax+lmax-1)*fac
+
+!--->    determine sin and cos of phi
+      s(0) = 0.0
+      s(1) = sph
+      c(0) = 1.0
+      c(1) = cph
+      cph2 = cph+cph
+      DO m = 2, lmax
+         s(m) = cph2*s(m-1)-s(m-2)
+         c(m) = cph2*c(m-1)-c(m-2)
+      ENDDO
+
+!--->    multiply in the normalization factors
+      DO l=0,lmax
+         ylm(l*(l+1)+1) = ynorm_dev(l*(l+1)+1)*cmplx(p(l,0),0.0)
+      ENDDO
+      DO m = 1, lmax
+         DO l = m, lmax
+            lm0 = l*(l+1)+1
+            ylms = p(l,m)*cmplx(c(m),s(m))
+            ylm(lm0+m) = ynorm_dev(lm0+m)*ylms
+            ylm(lm0-m) = conjg(ylms)*ynorm_dev(lm0-m)
+         ENDDO
+      ENDDO
+
+      END SUBROUTINE ylm4_dev
+#endif
+
+
       SUBROUTINE ylm4(lmax,v,ylm)
 !************************************************************
 !     generate the spherical harmonics for the vector v