LCOV - code coverage report
Current view: top level - math - ylm4.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 57 59 96.6 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             :       MODULE m_ylm
       2             :         IMPLICIT NONE
       3             :         PRIVATE
       4             : !---> check whether  or not normalizations are needed
       5             :         REAL, ALLOCATABLE, SAVE  :: ynorm(:)
       6             :         INTEGER,           SAVE  :: lmaxd = -1  ! initial value
       7             :         PUBLIC ylm4,ylmnorm_init
       8             : #ifdef CPP_GPU 
       9             :         PUBLIC ylm4_dev
      10             : #endif
      11             :       CONTAINS 
      12             : 
      13             : #ifdef CPP_GPU
      14             :       ATTRIBUTES(device) SUBROUTINE ylm4_dev(lmax,v,ylm)
      15             : !************************************************************
      16             : !     generate the spherical harmonics for the vector v
      17             : !     using a stable upward recursion in l.  (see notes
      18             : !     by m. weinert.)
      19             : !          m.weinert   january 1982
      20             : !     modified by R. Podloucky (added in ynorm); July 1989
      21             : !     cleaned up    mw 1995
      22             : !
      23             : !     modified to make use of f90 constructs. note that
      24             : !     the normalization is an internal subroutine and hence
      25             : !     can only be called from here. also, no need to dimension
      26             : !     arrays for ynorm, done dynamically.          mw 1999
      27             : !
      28             : !     GPU version added
      29             : !             U.Alekseeva      Oktober 2018
      30             : !************************************************************
      31             : !      USE m_juDFT
      32             :       INTEGER, VALUE, INTENT (IN) :: lmax
      33             :       REAL,    DEVICE, INTENT (IN) :: v(3)
      34             :       COMPLEX, DEVICE, INTENT (OUT):: ylm( (lmax+1)**2 )
      35             : 
      36             :       REAL, PARAMETER   :: small = 1.0e-12
      37             : 
      38             :       INTEGER  :: l,lm0,m
      39             :       REAL     :: fac,x,y,z,xy,r,rxy,cth,sth,cph,sph,cph2
      40             :       REAL,allocatable,device     :: c(:),s(:)
      41             :       REAL,allocatable,device     :: p(:,:)
      42             :       COMPLEX  :: ylms
      43             :       REAL,allocatable,device     :: ynorm_dev(:)
      44             :       REAL    a,cd,fpi
      45             : 
      46             :       allocate(c(0:lmax),s(0:lmax),ynorm_dev((lmax+1)**2))
      47             :       allocate(p(0:lmax,0:lmax))
      48             : 
      49             : !--->    calculate norm
      50             : 
      51             :       fpi = 4.0 * 3.1415926535897932  
      52             :       DO l=0,lmax
      53             :          lm0 = l*(l+1) + 1
      54             :          cd=1.0
      55             :          a = sqrt( (2*l+1)/fpi )
      56             :          ynorm_dev(lm0) = a
      57             :          DO m=1,l
      58             :             cd=cd/( (l+1-m)*(l+m) )
      59             :             ynorm_dev(lm0+m)  = a*sqrt(cd)
      60             :             ynorm_dev(lm0-m) = ( (-1.0)**m )*ynorm_dev(lm0+m)
      61             :          ENDDO
      62             :       ENDDO
      63             : 
      64             : !--->    calculate sin and cos of theta and phi
      65             :       x = v(1)
      66             :       y = v(2)
      67             :       z = v(3)
      68             :       xy  = x*x + y*y
      69             :       r   = sqrt(xy + z*z)
      70             :       rxy = sqrt(xy)
      71             : 
      72             :       IF (r.GT.small) THEN
      73             :          cth = z/r
      74             :          sth = rxy/r
      75             :       ELSE
      76             :          sth = 0.0
      77             :          cth = 1.0
      78             :       ENDIF
      79             :       IF(rxy.GT.small) THEN
      80             :          cph = x/rxy
      81             :          sph = y/rxy
      82             :       ELSE
      83             :          cph = 1.0
      84             :          sph = 0.0
      85             :       ENDIF
      86             : 
      87             : !---> generate associated legendre functions for m.ge.0
      88             :       fac = 1.0
      89             : !---> loop over m values
      90             :       DO m = 0, lmax - 1
      91             :          fac = -(m+m-1)*fac
      92             :          p(m,m) = fac
      93             :          p(m+1,m) = (m+m+1)*cth*fac
      94             : !--->    recurse upward in l
      95             :          DO l = m+2, lmax
      96             :             p(l,m)=((l+l-1)*cth*p(l-1,m)-(l+m-1)*p(l-2,m))/(l-m)
      97             :          ENDDO
      98             :          fac = fac*sth
      99             :       ENDDO
     100             :       p(lmax,lmax) = -(lmax+lmax-1)*fac
     101             : 
     102             : !--->    determine sin and cos of phi
     103             :       s(0) = 0.0
     104             :       s(1) = sph
     105             :       c(0) = 1.0
     106             :       c(1) = cph
     107             :       cph2 = cph+cph
     108             :       DO m = 2, lmax
     109             :          s(m) = cph2*s(m-1)-s(m-2)
     110             :          c(m) = cph2*c(m-1)-c(m-2)
     111             :       ENDDO
     112             : 
     113             : !--->    multiply in the normalization factors
     114             :       DO l=0,lmax
     115             :          ylm(l*(l+1)+1) = ynorm_dev(l*(l+1)+1)*cmplx(p(l,0),0.0)
     116             :       ENDDO
     117             :       DO m = 1, lmax
     118             :          DO l = m, lmax
     119             :             lm0 = l*(l+1)+1
     120             :             ylms = p(l,m)*cmplx(c(m),s(m))
     121             :             ylm(lm0+m) = ynorm_dev(lm0+m)*ylms
     122             :             ylm(lm0-m) = conjg(ylms)*ynorm_dev(lm0-m)
     123             :          ENDDO
     124             :       ENDDO
     125             : 
     126             :       END SUBROUTINE ylm4_dev
     127             : #endif
     128             : 
     129             : 
     130     6909592 :       SUBROUTINE ylm4(lmax,v,ylm)
     131             : !************************************************************
     132             : !     generate the spherical harmonics for the vector v
     133             : !     using a stable upward recursion in l.  (see notes
     134             : !     by m. weinert.)
     135             : !          m.weinert   january 1982
     136             : !     modified by R. Podloucky (added in ynorm); July 1989
     137             : !     cleaned up    mw 1995
     138             : !
     139             : !     modified to make use of f90 constructs. note that
     140             : !     the normalization is an internal subroutine and hence
     141             : !     can only be called from here. also, no need to dimension
     142             : !     arrays for ynorm, done dynamically.          mw 1999
     143             : !************************************************************
     144             :       USE m_juDFT
     145             :       INTEGER, INTENT (IN) :: lmax
     146             :       REAL,    INTENT (IN) :: v(3)
     147             :       COMPLEX, INTENT (OUT):: ylm( (lmax+1)**2 )
     148             : 
     149             :       REAL,    PARAMETER   :: small = 1.0e-12
     150             : 
     151             :       INTEGER l,lm0,m
     152             :       REAL    fac,x,y,z,xy,r,rxy,cth,sth,cph,sph,cph2
     153    13819184 :       REAL    p(0:lmax,0:lmax),c(0:max(1,lmax)),s(0:max(1,lmax))
     154             :       COMPLEX ylms
     155             : 
     156             : 
     157     6909592 :       IF (lmax.GT.lmaxd) THEN
     158           0 :          CALL juDFT_error("lmaxd too small in ylm4")
     159             :       ENDIF
     160             : 
     161             : !--->    calculate sin and cos of theta and phi
     162     6909592 :       x = v(1)
     163     6909592 :       y = v(2)
     164     6909592 :       z = v(3)
     165     6909592 :       xy  = x*x + y*y
     166     6909592 :       r   = sqrt(xy + z*z)
     167     6909592 :       rxy = sqrt(xy)
     168             : 
     169     6909592 :       IF (r.GT.small) THEN
     170     6909440 :          cth = z/r
     171     6909440 :          sth = rxy/r
     172             :       ELSE
     173             :          sth = 0.0
     174             :          cth = 1.0
     175             :       ENDIF
     176     6909592 :       IF(rxy.GT.small) THEN
     177     6805668 :          cph = x/rxy
     178     6805668 :          sph = y/rxy
     179             :       ELSE
     180             :          cph = 1.0
     181             :          sph = 0.0
     182             :       ENDIF
     183             : 
     184             : !---> generate associated legendre functions for m.ge.0
     185     6909592 :       fac = 1.0
     186             : !---> loop over m values
     187    63167056 :       DO m = 0, lmax - 1
     188    56257464 :          fac = -(m+m-1)*fac
     189    56257464 :          p(m,m) = fac
     190    56257464 :          p(m+1,m) = (m+m+1)*cth*fac
     191             : !--->    recurse upward in l
     192   267883780 :          DO l = m+2, lmax
     193   267883780 :             p(l,m)=((l+l-1)*cth*p(l-1,m)-(l+m-1)*p(l-2,m))/(l-m)
     194             :          ENDDO
     195    63167056 :          fac = fac*sth
     196             :       ENDDO
     197     6909592 :       p(lmax,lmax) = -(lmax+lmax-1)*fac
     198             : 
     199             : !--->    determine sin and cos of phi
     200     6909592 :       s(0) = 0.0
     201     6909592 :       s(1) = sph
     202     6909592 :       c(0) = 1.0
     203     6909592 :       c(1) = cph
     204     6909592 :       cph2 = cph+cph
     205    56257464 :       DO m = 2, lmax
     206    49347872 :          s(m) = cph2*s(m-1)-s(m-2)
     207    56257464 :          c(m) = cph2*c(m-1)-c(m-2)
     208             :       ENDDO
     209             : 
     210             : !--->    multiply in the normalization factors
     211   133243704 :       DO l=0,lmax
     212    70076648 :          ylm(l*(l+1)+1) = ynorm(l*(l+1)+1)*cmplx(p(l,0),0.0)
     213             :       ENDDO
     214   119424520 :       DO m = 1, lmax
     215   598934616 :          DO l = m, lmax
     216   267883780 :             lm0 = l*(l+1)+1
     217   267883780 :             ylms = p(l,m)*cmplx(c(m),s(m))
     218   267883780 :             ylm(lm0+m) = ynorm(lm0+m)*ylms
     219   324141244 :             ylm(lm0-m) = conjg(ylms)*ynorm(lm0-m)
     220             :          ENDDO
     221             :       ENDDO
     222             : 
     223     6909592 :       RETURN
     224             : 
     225             :       END SUBROUTINE ylm4
     226             : 
     227         128 :       SUBROUTINE ylmnorm_init(lmax)
     228             : !********************************************************************
     229             : !     normalization constants for ylm (internal subroutine has access
     230             : !     to lmax and ynorm from above)
     231             : !********************************************************************
     232             : !$    USE OMP_LIB
     233             :       use m_juDFT
     234             :       USE m_constants, ONLY : pimach
     235             :       INTEGER,INTENT(IN):: lmax
     236             :       INTEGER l,lm0,m
     237             :       REAL    a,cd,fpi
     238         128 :       IF ( allocated(ynorm) ) DEALLOCATE(ynorm)
     239         128 :       ALLOCATE ( ynorm( (lmax+1)**2 ) )   ! allocate array
     240         128 :       lmaxd = lmax
     241             : 
     242         128 :       fpi = 4.0*pimach()
     243         128 : !$    if (omp_in_parallel()) call juDFT_error(  &
     244           0 : !$          "ylmnorm should not called in parallel",calledby="ylm4.f")
     245        2064 :       DO l=0,lmax
     246        1936 :          lm0 = l*(l+1) + 1
     247        1936 :          cd=1.0
     248        1936 :          a = sqrt( (2*l+1)/fpi )
     249        1936 :          ynorm(lm0) = a
     250       20520 :          DO m=1,l
     251       18584 :             cd=cd/( (l+1-m)*(l+m) )
     252       18584 :             ynorm(lm0+m)  = a*sqrt(cd)
     253       20520 :             ynorm(lm0-m) = ( (-1.0)**m )*ynorm(lm0+m)
     254             :          ENDDO
     255             :       ENDDO
     256             : 
     257         128 :       RETURN
     258             :       END SUBROUTINE ylmnorm_init
     259             : 
     260             :       END MODULE m_ylm

Generated by: LCOV version 1.13