LCOV - code coverage report
Current view: top level - math - intgr.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 124 146 84.9 %
Date: 2019-09-08 04:53:50 Functions: 6 9 66.7 %

          Line data    Source code
       1             : MODULE m_intgr
       2             : 
       3             :   !**********************************************************************
       4             :   ! intgr[0-3]:
       5             :   !    Integrators of a function y(jri) on a logarithmic mesh with jri 
       6             :   !    mesh points. The output is either a scalar [z] or a field [z(jri)].
       7             :   !    Either the first meshpoint [r0,r1] and the increment [h] or the 
       8             :   !    array of meshpoints [rmsh(jri)] is supplied:
       9             :   !
      10             :   ! intgz0 & intgz1 :
      11             :   !     (in)definite integrators on linear mesh. tail=.true. will include
      12             :   !     a tail correction assuming that the function is a simple
      13             :   !     decaying exponential between the first mesh point and infinity.
      14             :   !     y contains the nmz function values tabulated at a spacing of h.
      15             :   !
      16             :   !            integrator:      ---- input ----      output
      17             :   !    intgr0:   definite       y r0        h jri  | z
      18             :   !    intgr1: indefinite       y r1        h jri  | z(jri)
      19             :   !    intgr2: indefinite       y rmsh(jri) h jri  | z(jri)
      20             :   !    intgr3:   definite       y rmsh(jri) h jri  | z
      21             :   !    intgz0:   definite       y tail      n nmz  ! z
      22             :   !    intgz1: indefinite       y tail      n nmz  ! z(nmz)
      23             :   !
      24             :   !                                                            m. weinert
      25             :   !**********************************************************************
      26             : 
      27             : #include"cpp_double.h"
      28             : 
      29             :   IMPLICIT NONE
      30             : 
      31             :   !INTRINSIC exp,log
      32             :   INTERFACE 
      33             :     REAL FUNCTION CPP_BLAS_sdot( n, f1, is1, f2, is2 )
      34             :       INTEGER, INTENT (IN) :: n, is1, is2
      35             :       REAL,    INTENT (IN) :: f1(n), f2(n)
      36             :     END FUNCTION 
      37             :   END INTERFACE
      38             : 
      39             : !  interface intgz1
      40             : !    module procedure intgz1Real, intgz1Complex
      41             : !  end interface
      42             : 
      43             :   interface intgz1Reverse
      44             :     module procedure intgz1RealReverse, intgz1ComplexReverse
      45             :   end interface
      46             :   
      47             : 
      48             :   INTEGER, PARAMETER, PRIVATE :: nr  = 7 , nr1 = 6
      49             :   REAL,    PARAMETER, PRIVATE :: h0 = 140. , zero = 0.0e0
      50             :   
      51             :   ! lagrangian integration coefficients (simpson 7 point rule: error  h**9)
      52             :   
      53             :   INTEGER, DIMENSION(7),   PARAMETER, PRIVATE :: ih = (/41,216,27,272,27,216,41/)
      54             : 
      55             :   REAL,    DIMENSION(7,5), PARAMETER, PRIVATE :: a = RESHAPE( &
      56             :        (/19087.,65112.,-46461., 37504.,-20211., 6312.,-863., &
      57             :           -863.,25128., 46989.,-16256.,  7299.,-2088., 271., &
      58             :            271.,-2760., 30819., 37504., -6771., 1608.,-191., &
      59             :           -191., 1608., -6771., 37504., 30819.,-2760., 271., &
      60             :            271.,-2088.,  7299.,-16256., 46989.,25128.,-863./),(/7,5/))
      61             : 
      62             : 
      63             :   CONTAINS
      64             : 
      65             :   !**********************************************************************
      66     2222395 :   SUBROUTINE intgr0( y, r0, h, jri, z )
      67             :     !**********************************************************************
      68             :     !     ..
      69             :     !     .. Arguments ..
      70             :     INTEGER,    INTENT (IN)  :: jri
      71             :     REAL,       INTENT (IN)  :: h, r0
      72             :     REAL,       INTENT (IN)  :: y(jri)
      73             :     REAL,       INTENT (OUT) :: z
      74             :     !     ..
      75             :     !     .. Locals ..
      76             :     INTEGER :: m, n0, nsteps
      77             :     REAL    :: dr, r(7)
      78             :     INTEGER :: i, j
      79             :     REAL    :: alpha, yr(7)
      80             : 
      81             :     !
      82             :     !--->    integral from 0 to r1 approximated by leading term in power
      83             :     !--->    series expansion of y(r)
      84             :     !
      85     2222395 :     z = zero
      86     2222395 :     IF (y(1)*y(2).GT.zero) THEN
      87     2192465 :       alpha = 1.0 + log(y(2)/y(1))/h
      88     2192465 :       IF (alpha.GT.zero) z = r0*y(1)/alpha
      89             :     ENDIF
      90             :     !
      91             :     !--->    determine steps and starting point for simpson
      92             :     !
      93     2222395 :     nsteps = (jri-1)/nr1
      94     2222395 :     n0 = jri - nr1*nsteps
      95     2222395 :     dr = exp(h)
      96     2222395 :     r(1) = r0
      97    15556765 :     DO i = 2,7
      98    15556765 :       r(i) = dr*r(i-1)
      99             :     ENDDO 
     100             :     !
     101             :     !--->    lagrange integration for points 1<j<n0, error: h**9
     102             :     !
     103     2222395 :     IF (n0.GT.1) THEN
     104     5155005 :       DO i = 1,7
     105     2749336 :         yr(i) = r(i)*y(i)
     106             :       ENDDO
     107     2699655 :       DO j = 1,n0 - 1
     108     1521661 :         z = z + h*CPP_BLAS_sdot(7,a(1,j),1,yr,1)/60480.
     109             :       ENDDO
     110             :     ENDIF
     111     2222395 :     r(1) = r(n0)
     112             :     !
     113             :     !--->    simpson integration
     114             :     !
     115   231934790 :     DO m = 1,nsteps
     116  2986261135 :       DO i = 2,nr
     117  1607986765 :         r(i) = dr*r(i-1)
     118             :       ENDDO
     119  3445685925 :       DO i = 1,nr
     120  1837699160 :         yr(i) = h*ih(i)*r(i)/h0
     121             :       ENDDO
     122   229712395 :       z = z + CPP_BLAS_sdot(nr,yr,1,y(n0),1)
     123   229712395 :       n0 = n0 + nr1
     124   231934790 :       r(1) = r(nr)
     125             :     ENDDO
     126             : 
     127     2222395 :     RETURN
     128             :   END SUBROUTINE intgr0
     129             : 
     130             :   !**********************************************************************
     131       16999 :   SUBROUTINE intgr1( y, r1, h, jri, z )
     132             :     !**********************************************************************
     133             :     !     .. 
     134             :     !     .. Arguments ..
     135             :     INTEGER, INTENT (IN) :: jri
     136             :     REAL,    INTENT (IN) :: h,r1
     137             :     REAL,    INTENT (IN) :: y(jri)
     138             :     REAL,    INTENT (OUT):: z(jri)
     139             :     !     ..
     140             :     !     .. Locals ..
     141             :     REAL    :: dr, rr, r(7)
     142             :     INTEGER :: i, j
     143             :     REAL    :: alpha, yr(7)
     144             :     !
     145             :     !--->    integral from 0 to r1 approximated by leading term in power
     146             :     !--->    series expansion of y(r)
     147             :     !
     148       16999 :     z(1) = zero
     149       16999 :     IF (y(1)*y(2).GT.zero) THEN
     150       16999 :       alpha = 1.0 + log(y(2)/y(1))/h
     151       16999 :       IF (alpha.GT.zero) z(1) = r1*y(1)/alpha
     152             :     ENDIF
     153             :     !
     154             :     !--->    lagrange integration for points 1<j<nr, error: h**9
     155             :     !
     156       16999 :     dr = exp(h)
     157       16999 :     rr = r1
     158      135992 :     DO i = 1,7
     159      118993 :       yr(i) = rr*y(i)
     160      118993 :       r(i) = rr
     161      135992 :       rr = dr*rr
     162             :     ENDDO
     163      186989 :     DO j = 1,nr - 2
     164      101994 :       z(j+1) = z(j) + h*CPP_BLAS_sdot(7,a(1,j),1,yr,1)/60480.
     165             :     ENDDO
     166             :     !
     167             :     !--->    simpson integration, j>nr-1
     168             :     !
     169      254985 :     DO i = 1,nr
     170      135992 :       r(i) = h*ih(i)*r(i)/h0
     171             :     ENDDO
     172    27703015 :     DO j = nr,jri
     173    13843008 :       z(j) = z(j-nr1) + CPP_BLAS_sdot(nr,r,1,y(j-nr1),1)
     174   110761063 :       DO i = 1,7
     175   110744064 :         r(i) = dr*r(i)
     176             :       ENDDO
     177             :     ENDDO
     178             : 
     179       16999 :     RETURN
     180             :   END SUBROUTINE intgr1
     181             : 
     182             :   !**********************************************************************
     183       11046 :   SUBROUTINE intgr2( y, rmsh, h, jri, z )
     184             :     !**********************************************************************
     185             :     !     ..
     186             :     !     .. Arguments ..
     187             :     INTEGER, INTENT (IN) :: jri
     188             :     REAL,    INTENT (IN) :: h
     189             :     REAL,    INTENT (IN) :: rmsh(jri), y(jri)
     190             :     REAL,    INTENT (OUT):: z(jri)
     191             :     !     ..
     192             :     !     .. Locals ..
     193             :     REAL    :: dr, r(7)
     194             :     INTEGER :: i, j
     195             :     REAL    :: alpha, yr(7)
     196             :     !
     197             :     !--->    integral from 0 to r1 approximated by leading term in power
     198             :     !--->    series expansion of y(r)
     199             :     !
     200       11046 :     z(1) = zero
     201       11046 :     IF (y(1)*y(2).GT.zero) THEN 
     202       10292 :       alpha = 1.0 + log(y(2)/y(1))/h
     203       10292 :       IF (alpha.GT.zero) z(1) = rmsh(1)*y(1)/alpha
     204             :     ENDIF
     205             :     !
     206             :     !--->    lagrange integration for points 1<j<nr, error: h**9
     207             :     !
     208       11046 :     dr = exp(h)
     209       88368 :     DO i = 1,7
     210       77322 :       r(i) = rmsh(i)
     211       88368 :       yr(i) = rmsh(i)*y(i)
     212             :     ENDDO
     213      121506 :     DO j = 1,nr - 2
     214       66276 :       z(j+1) = z(j) + h*CPP_BLAS_sdot(7,a(1,j),1,yr,1)/60480.
     215             :     ENDDO
     216             :     !
     217             :     !--->    simpson integration, j>nr-1
     218             :     !
     219      165690 :     DO i = 1,nr
     220       88368 :       r(i) = h*ih(i)*r(i)/h0
     221             :     ENDDO
     222    13584794 :     DO j = nr,jri
     223     6786874 :       z(j) = z(j-nr1) + CPP_BLAS_sdot(nr,r,1,y(j-nr1),1)
     224    54306038 :       DO i = 1,7
     225    54294992 :         r(i) = dr*r(i)
     226             :       ENDDO
     227             :     ENDDO
     228             : 
     229       11046 :     RETURN
     230             :   END SUBROUTINE intgr2
     231             : 
     232             :   !**********************************************************************
     233     1212843 :   SUBROUTINE intgr3( y, r, h, jri, z )
     234             :     !**********************************************************************
     235             :     !     ..
     236             :     !     .. Arguments ..
     237             :     INTEGER, INTENT (IN) :: jri
     238             :     REAL,    INTENT (IN) :: h
     239             :     REAL,    INTENT (IN) :: r(jri)
     240             :     REAL,    INTENT (IN) :: y(jri)
     241             :     REAL,    INTENT (OUT):: z
     242             :     !     ..
     243             :     !     .. Locals ..
     244             :     INTEGER :: m, n0, nsteps
     245             :     REAL    :: tiny, yr(nr), h1, z1, ih1(nr)
     246             :     INTEGER :: i, j
     247             :     REAL    :: alpha
     248             :     !
     249             :     !--->    integral from 0 to r1 approximated by leading term in power
     250             :     !--->    series expansion of y(r)
     251             :     !
     252             :     !      DO i=1,jri
     253             :     !        IF (abs(y(i)).LT.tiny) y(i) = tiny
     254             :     !      ENDDO
     255             :     !
     256     1212843 :     z = zero
     257     1212843 :     IF (y(1)*y(2).GT.zero) THEN
     258      937845 :       alpha = 1.0 + log(y(2)/y(1))/h
     259      937845 :       IF (alpha.GT.zero) z = r(1)*y(1)/alpha
     260             :     ENDIF
     261             :     !
     262             :     !--->    determine steps and starting point for simpson
     263             :     !
     264     1212843 :     nsteps = (jri-1)/nr1
     265     1212843 :     n0 = jri - nr1*nsteps
     266             :     !
     267             :     !--->    lagrange integration for points 1<j<n0, error: h**9
     268             :     !
     269     1212843 :     IF (n0.GT.1) THEN
     270     3843150 :       DO i = 1,7
     271     2049680 :         yr(i) = r(i)*y(i)
     272             :       ENDDO
     273             :       z1 = 0.
     274     2288870 :       DO j = 1,n0 - 1
     275     1272540 :         z1 = z1 + CPP_BLAS_sdot(7,a(1,j),1,yr,1)
     276             :       ENDDO
     277      256210 :       z = z + z1 * h / 60480.
     278             :     END IF
     279             :     !
     280             :     !--->    simpson integration
     281             :     !
     282     1212843 :     h1 = h / h0
     283     9702744 :     DO i = 1,nr
     284     9702744 :       ih1(i) = h1 * ih(i)
     285             :     ENDDO
     286   248980529 :     DO m = 1,nsteps
     287  1858257645 :       DO i = 1,nr
     288   991070744 :         yr(i) = ih1(i)*r(i+n0-1)
     289             :       ENDDO
     290   123883843 :       z = z + CPP_BLAS_sdot(nr,yr,1,y(n0),1)
     291   125096686 :       n0 = n0 + nr1
     292             :     ENDDO
     293             : 
     294     1212843 :     RETURN
     295             :   END SUBROUTINE intgr3
     296             : 
     297             :   !**********************************************************************
     298     1196231 :   SUBROUTINE intgz0( y, h, nmz, z, tail )
     299             :     !**********************************************************************
     300             :     !     ..
     301             :     !     .. Arguments ..
     302             :     INTEGER, INTENT (IN)  :: nmz
     303             :     REAL,    INTENT (IN)  :: h
     304             :     REAL,    INTENT (IN)  :: y(nmz)
     305             :     LOGICAL, INTENT (IN)  :: tail
     306             :     REAL,    INTENT (OUT) :: z
     307             :     !     ..
     308             :     !     .. Locals ..
     309             :     REAL    :: yl, ys
     310             :     INTEGER :: m, n0, nsteps
     311             :     INTEGER :: i, j
     312             :     REAL    :: alpha, yr(7)
     313             :     !
     314             :     !--->    integral from minus infinity to the first mesh point assuming
     315             :     !--->    exponential decay in this region. this contribution is limited
     316             :     !--->    to alpha>0.1 which corresponds to square of wavefunctions
     317             :     !--->    of energy > 0.00125 a.u. below the vacuum level
     318             :     !
     319     1196231 :     z = zero
     320     1196231 :     IF (tail) THEN
     321      864732 :       IF (y(1)*y(2).GT.zero) THEN
     322      755178 :         alpha = log(y(2)/y(1))/h
     323      755178 :         IF (alpha.GT.0.1) z = y(1)/alpha
     324             :       ENDIF
     325             :     ENDIF
     326             :     !
     327             :     !--->    determine steps and starting point for simpson
     328             :     !
     329     1196231 :     nsteps = (nmz-1)/nr1
     330     1196231 :     n0 = nmz - nr1*nsteps
     331             :     !
     332             :     !--->    lagrange integration for points 1<j<n0, error: h**9
     333             :     !
     334     1196231 :     yl = zero
     335     1196231 :     IF (n0.GT.1) THEN
     336     7867052 :       DO j = 1, n0 - 1
     337     4508066 :         yl = yl + CPP_BLAS_sdot(7,a(1,j),1,y,1)
     338             :       ENDDO
     339     1149080 :       yl = h*yl/60480.
     340             :     END IF
     341             :     !
     342             :     !--->    simpson integration
     343             :     !
     344     1196231 :     ys = zero
     345     1196231 :     n0 = n0 - 1
     346    19042157 :     DO m = 1,nsteps
     347   267688890 :       DO i = 1,nr
     348   142767408 :         ys = ys + ih(i)*y(n0+i)
     349             :       ENDDO
     350    19042157 :       n0 = n0 + nr1
     351             :     ENDDO
     352     1196231 :     ys = h*ys/h0
     353     1196231 :     z = z + yl + ys
     354             : 
     355     1196231 :     RETURN
     356             :   END SUBROUTINE intgz0
     357             : 
     358             :   !**********************************************************************
     359        2240 :   SUBROUTINE intgz1( y, h, nmz, z, tail )
     360             :     !**********************************************************************
     361             :     !     .. 
     362             :     !     .. Arguments ..
     363             :     INTEGER, INTENT (IN)  :: nmz
     364             :     LOGICAL, INTENT (IN)  :: tail
     365             :     REAL,    INTENT (IN)  :: h
     366             :     REAL,    INTENT (IN)  :: y(nmz)
     367             :     REAL,    INTENT (OUT) :: z(nmz)
     368             :     !     ..
     369             :     !     .. Locals ..
     370             :     REAL    :: yl, ys
     371             :     REAL, PARAMETER :: eps = 1.e-38
     372             :     INTEGER :: i, j
     373             :     REAL    :: alpha, yr(7)
     374             :     !
     375             :     !--->    integral from minus infinity to the first mesh point assuming
     376             :     !--->    exponential decay in this region. this contribution is limited
     377             :     !--->    to alpha>0.1 which corresponds to square of wavefunctions
     378             :     !--->    of energy > 0.00125 a.u. below the vacuum level
     379             :     !
     380        2240 :     z(1) = zero
     381        2240 :     IF (tail) THEN
     382        2240 :       IF (abs(y(1)).GT.eps) THEN
     383          28 :         IF (y(1)*y(2).GT.zero) THEN
     384          28 :           alpha = log(y(2)/y(1))/h
     385          28 :           IF (alpha.GT.0.1) z(1) = y(1)/alpha
     386             :         ENDIF
     387             :       ENDIF
     388             :     ENDIF
     389             :     !
     390             :     !--->    lagrange integration for points 1<j<n0, error: h**9
     391             :     !
     392       24640 :     DO j = 1,nr - 2
     393       11200 :       yl = 0
     394       11200 :       yl = yl + CPP_BLAS_sdot(7,a(1,j),1,y,1)
     395       13440 :       z(j+1) = z(j) + h*yl/60480.
     396             :     ENDDO
     397             :     !
     398             :     !--->    simpson integration
     399             :     !
     400      423360 :     DO j = nr,nmz
     401             :       ys = zero
     402     3158400 :       DO i = 1,nr
     403     1684480 :         ys = ys + ih(i)*y(j-nr+i)
     404             :       ENDDO
     405      212800 :       z(j) = z(j-nr1) + h*ys/h0
     406             :     ENDDO
     407             : 
     408        2240 :     RETURN
     409             :   END SUBROUTINE intgz1
     410             : 
     411             : 
     412             : 
     413           0 :   subroutine intgz1Complex( y, h, nmz, z, tail )
     414             : 
     415             :     integer, intent(in)  :: nmz
     416             :     logical, intent(in)  :: tail
     417             :     real,    intent(in)  :: h
     418             :     complex, intent(in)  :: y(nmz)
     419             :     complex, intent(out) :: z(nmz)
     420             : 
     421           0 :     real                 :: zr(nmz), zi(nmz)
     422             : 
     423           0 :     call intgz1(  real(y), h, nmz, zr, tail )
     424           0 :     call intgz1( aimag(y), h, nmz, zi, tail )
     425           0 :     z = cmplx( zr, zi )
     426             : 
     427           0 :   end subroutine intgz1Complex
     428             : 
     429             : 
     430             : 
     431           0 :   subroutine intgz1RealReverse( y, h, nmz, z, tail )
     432             : 
     433             :     integer, intent(in)  :: nmz
     434             :     logical, intent(in)  :: tail
     435             :     real,    intent(in)  :: h
     436             :     real,    intent(in)  :: y(nmz)
     437             :     real,    intent(out) :: z(nmz)
     438             : 
     439           0 :     real                 :: y_reverse(nmz), z_reverse(nmz)
     440             :     integer              :: i
     441             : 
     442           0 :     do i = 1, nmz
     443           0 :       y_reverse(i) = y(nmz+1-i)
     444             :     end do
     445           0 :     call intgz1( y_reverse, h, nmz, z_reverse, tail )
     446           0 :     do i = 1, nmz
     447           0 :       z(i) = z_reverse(nmz+1-i)
     448             :     end do    
     449             : 
     450           0 :   end subroutine intgz1RealReverse
     451             : 
     452             : 
     453             : 
     454           0 :   subroutine intgz1ComplexReverse( y, h, nmz, z, tail )
     455             : 
     456             :     integer, intent(in)  :: nmz
     457             :     logical, intent(in)  :: tail
     458             :     real,    intent(in)  :: h
     459             :     complex, intent(in)  :: y(nmz)
     460             :     complex, intent(out) :: z(nmz)
     461             : 
     462           0 :     complex              :: y_reverse(nmz), z_reverse(nmz)
     463             :     integer              :: i
     464             : 
     465           0 :     do i = 1, nmz
     466           0 :       y_reverse(i) = y(nmz+1-i)
     467             :     end do
     468           0 :     call intgz1Complex( y_reverse, h, nmz, z_reverse, tail )
     469           0 :     do i = 1, nmz
     470           0 :       z(i) = z_reverse(nmz+1-i)
     471             :     end do    
     472             : 
     473           0 :   end subroutine intgz1ComplexReverse
     474             : 
     475             : 
     476             : 
     477             : END MODULE m_intgr

Generated by: LCOV version 1.13