LCOV - code coverage report
Current view: top level - math - grule.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 30 30 100.0 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_grule
       2             :       CONTAINS 
       3         514 :       SUBROUTINE grule(n,x,w)
       4             : c***********************************************************************
       5             : c     determines the (n+1)/2 nonnegative points x(i) and
       6             : c     the corresponding weights w(i) of the n-point
       7             : c     gauss-legendre integration rule, normalized to the
       8             : c     interval (-1,1). the x(i) appear in descending order.
       9             : c     this routine is from 'methods of numerical integration',
      10             : c     p.j. davis and p. rabinowitz, page 369.
      11             : c                                                            m.w.
      12             : c***********************************************************************
      13             : 
      14             :       USE m_constants
      15             :       IMPLICIT NONE
      16             : C     ..
      17             : C     .. Arguments ..
      18             :       INTEGER, INTENT (IN) :: n
      19             :       REAL,    INTENT(OUT) :: w(n/2),x(n/2)
      20             : C     ..
      21             : C     .. Locals ..
      22             :       INTEGER i,it,k,m
      23             :       REAL d1,d2pn,d3pn,d4pn,den,dp,dpn,e1,fx,h
      24             :       REAL p,pk,pkm1,pkp1,t,t1,u,v,x0
      25             : C     ..
      26             : C     ..
      27         514 :       m = (n+1)/2
      28         514 :       e1 = n* (n+1)
      29             : 
      30        3404 :       DO i = 1,m
      31        2890 :          t = (4*i-1)*pi_const/ (4*n+2)
      32        2890 :          x0 = (1.- (1.-1./n)/ (8.*n*n))*cos(t)
      33             : c--->    iterate on the value  (m.w. jan. 1982)
      34        8670 :          DO it = 1,2
      35             :             pkm1 = 1.
      36             :             pk = x0
      37      131388 :             DO k = 2,n
      38       62804 :                t1 = x0*pk
      39       62804 :                pkp1 = t1 - pkm1 - (t1-pkm1)/k + t1
      40       62804 :                pkm1 = pk
      41       68584 :                pk = pkp1
      42             :             ENDDO
      43        5780 :             den = 1. - x0*x0
      44        5780 :             d1 = n* (pkm1-x0*pk)
      45        5780 :             dpn = d1/den
      46        5780 :             d2pn = (2.*x0*dpn-e1*pk)/den
      47        5780 :             d3pn = (4.*x0*d2pn+ (2.-e1)*dpn)/den
      48        5780 :             d4pn = (6.*x0*d3pn+ (6.-e1)*d2pn)/den
      49        5780 :             u = pk/dpn
      50        5780 :             v = d2pn/dpn
      51        5780 :             h = -u* (1.+.5*u* (v+u* (v*v-u*d3pn/ (3.*dpn))))
      52        5780 :             p = pk + h* (dpn+.5*h* (d2pn+h/3.* (d3pn+.25*h*d4pn)))
      53        5780 :             dp = dpn + h* (d2pn+.5*h* (d3pn+h*d4pn/3.))
      54        5780 :             h = h - p/dp
      55        8670 :             x0 = x0 + h
      56             :          ENDDO
      57        2890 :          x(i) = x0
      58             :          fx = d1 - h*e1* (pk+.5*h* (dpn+h/3.* (d2pn+.25*h* (d3pn+
      59        2890 :      +        .2*h*d4pn))))
      60        3404 :          w(i) = 2.* (1.-x(i)*x(i))/ (fx*fx)
      61             :       ENDDO
      62             : 
      63         514 :       IF (m+m.GT.n) x(m) = 0.
      64         514 :       RETURN
      65             :       END SUBROUTINE grule
      66             :       END MODULE M_grule

Generated by: LCOV version 1.13