- Added preliminary implementation for Legendre polynomials

- Added preliminary implementation for Legendre polynomial first derivatives
- Added preliminary implementation for Gauss-Legendre quadrature points/weights
- Added preliminary implementation for Gauss-Legendre-Lobatto quadrature points/weights

Author: hsnyder

src/CMakeLists.txt

@@ -41,6 +41,9 @@ set(SRC
     stdlib_kinds.f90
     stdlib_logger.f90
     stdlib_system.F90
+    stdlib_functions.f90
+    stdlib_functions_legendre.f90
+    stdlib_quadrature_gauss.f90
     ${outFiles}
 )
 

src/Makefile.manual

@@ -4,6 +4,8 @@ SRC = f18estop.f90 \
       stdlib_bitsets_64.f90 \
       stdlib_bitsets_large.f90 \
       stdlib_error.f90 \
+      stdlib_functions.f90 \
+      stdlib_functions_legendre.f90 \
       stdlib_io.f90 \
       stdlib_kinds.f90 \
       stdlib_linalg.f90 \
@@ -12,6 +14,7 @@ SRC = f18estop.f90 \
       stdlib_optval.f90 \
       stdlib_quadrature.f90 \
       stdlib_quadrature_trapz.f90 \
+      stdlib_quadrature_gauss.f90 \
       stdlib_stats.f90 \
       stdlib_stats_mean.f90 \
       stdlib_stats_moment.f90 \
@@ -50,6 +53,8 @@ stdlib_bitsets.o: stdlib_kinds.o
 stdlib_bitsets_64.o: stdlib_bitsets.o
 stdlib_bitsets_large.o: stdlib_bitsets.o
 stdlib_error.o: stdlib_optval.o
+stdlib_functions.o: stdlib_kinds.o
+stdlib_functions_legendre.o: stdlib_kinds.o
 stdlib_io.o: \
 	stdlib_error.o \
 	stdlib_optval.o \
@@ -58,6 +63,7 @@ stdlib_linalg_diag.o: stdlib_kinds.o
 stdlib_logger.o: stdlib_ascii.o stdlib_optval.o
 stdlib_optval.o: stdlib_kinds.o
 stdlib_quadrature.o: stdlib_kinds.o
+stdlib_quadrature_gauss.o: stdlib_kinds.o
 stdlib_stats_mean.o: \
 	stdlib_optval.o \
 	stdlib_kinds.o \

src/stdlib_functions.f90

@@ -0,0 +1,34 @@
+module stdlib_functions
+    use stdlib_kinds, only: sp, dp, qp
+
+    implicit none
+
+    private
+
+    public :: legendre 
+    public :: dlegendre 
+
+
+    interface legendre
+        !! version: experimental
+        !! 
+        !! Legendre polynomial
+        pure elemental module function legendre_fp64(N,x) result(L)
+            integer, intent(in) :: N
+            real(dp), intent(in) :: x
+            real(dp) :: L
+        end function
+    end interface
+
+    interface dlegendre
+        !! version: experimental
+        !! 
+        !! First derivative Legendre polynomial
+        pure elemental module function dlegendre_fp64(N,x) result(DL)
+            integer, intent(in) :: N
+            real(dp), intent(in) :: x
+            real(dp) :: DL
+        end function
+    end interface
+
+end module stdlib_functions

src/stdlib_functions_legendre.f90

@@ -0,0 +1,73 @@
+submodule (stdlib_functions) stdlib_functions_legendre
+    use stdlib_kinds, only: sp, dp, qp
+    implicit none
+
+contains
+
+    ! derivatives of Legendre polynomials
+    ! unspecified behaviour if N is negative
+    pure elemental module function dlegendre_fp64(N,x) result(DL)
+        integer, intent(in) :: N
+        real(dp), intent(in) :: x
+        real(dp) :: DL
+
+        select case(N)
+            case(0)
+                DL = 0
+            case(1)
+                DL = 1
+            case default
+                block
+                    real(dp) :: L_down1, L_down2, L
+                    real(dp) :: DL_down1, DL_down2
+                    integer :: i 
+
+                    L_down1  = x
+                    DL_down1 = 1
+
+                    L_down2  = 1
+                    DL_down2 = 0
+
+                    do i = 2, N
+                        L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
+                        DL = DL_down2 + (2*i-1)*L_down1
+                        L_down2 = L_down1
+                        L_down1 = L
+                        DL_down2 = DL_down1
+                        DL_down1 = DL
+                    end do
+                end block
+        end select
+    end function
+
+    ! Legendre polynomials
+    ! unspecified behaviour if N is negative
+    pure elemental module function legendre_fp64(N,x) result(L)
+        integer, intent(in) :: N
+        real(dp), intent(in) :: x
+        real(dp) :: L
+        select case(N)
+            case(0)
+                L  = 1
+            case(1)
+                L  = x
+            case default
+                block
+                    real(dp) :: L_down1, L_down2
+                    integer :: i 
+
+                    L_down1  = x
+                    L_down2  = 1
+
+                    do i = 2, N
+                        L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
+                        L_down2 = L_down1
+                        L_down1 = L
+                    end do
+                end block
+        end select
+    end function
+
+end submodule
+
+

src/stdlib_quadrature.fypp

@@ -12,6 +12,8 @@ module stdlib_quadrature
     public :: trapz_weights
     public :: simps
     public :: simps_weights
+    public :: gauss_legendre
+    public :: gauss_legendre_lobatto
 
 
     interface trapz
@@ -90,6 +92,28 @@ module stdlib_quadrature
     end interface simps_weights
 
 
+    interface gauss_legendre
+        !! version: experimental
+        !!
+        !! Computes Gauss-Legendre quadrature nodes and weights.
+        pure module subroutine gauss_legendre_fp64 (x, w, interval)
+            real(dp), intent(out) :: x(:), w(:)
+            real(dp), intent(in), optional :: interval(2)
+        end subroutine
+    end interface gauss_legendre
+
+
+    interface gauss_legendre_lobatto
+        !! version: experimental
+        !!
+        !! Computes Gauss-Legendre-Lobatto quadrature nodes and weights.
+        pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
+            real(dp), intent(out) :: x(:), w(:)
+            real(dp), intent(in), optional :: interval(2)
+        end subroutine
+    end interface gauss_legendre_lobatto
+
+
     ! Interface for a simple f(x)-style integrand function.
     ! Could become fancier as we learn about the performance
     ! ramifications of different ways to do callbacks.
@@ -103,4 +127,6 @@ module stdlib_quadrature
         #:endfor
     end interface
 
+
+
 end module stdlib_quadrature

src/stdlib_quadrature_gauss.f90

@@ -0,0 +1,115 @@
+submodule (stdlib_quadrature) stdlib_quadrature_gauss
+    use stdlib_kinds, only: sp, dp, qp
+    use stdlib_functions, only: legendre, dlegendre
+    implicit none
+
+    real(dp), parameter :: PI = acos(-1._dp)
+    real(dp), parameter :: tolerance = 4._dp * epsilon(1._dp)
+    integer, parameter :: newton_iters = 100
+
+contains
+
+    pure module subroutine gauss_legendre_fp64 (x, w, interval)
+        real(dp), intent(out) :: x(:), w(:)
+        real(dp), intent(in), optional :: interval(2)
+
+        associate (N => size(x)-1 )
+        select case (N)
+            case (0)
+                x = 0._dp
+                w = 2._dp
+            case (1)
+                x = [-sqrt(1._dp/3._dp), sqrt(1._dp/3._dp)]
+                w = [1._dp, 1._dp]
+            case default
+                block
+                integer :: i,j
+                real(dp) :: leg, dleg, delta
+
+                do i = 0, int(floor((N+1)/2._dp)-1)
+                    x(i+1) = -cos((2*i+1)/(2._dp*N+2._dp) * PI)
+                    do j = 0, newton_iters-1
+                        leg  = legendre(N+1,x(i+1))
+                        dleg = dlegendre(N+1,x(i+1))
+                        delta = -leg/dleg
+                        x(i+1) = x(i+1) + delta
+                        if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
+                    end do
+                    x(N-i+1) = -x(i+1)
+
+                    dleg = dlegendre(N+1,x(i+1))
+                    w(i+1)   = 2._dp/((1-x(i+1)**2)*dleg**2) 
+                    w(N-i+1) = w(i+1)
+                end do
+
+                if (mod(N,2) == 0) then
+                    x(N/2+1) = 0.0
+
+                    dleg = dlegendre(N+1, 0.0_dp)
+                    w(N/2+1) = 2._dp/(dleg**2) 
+                end if
+                end block
+        end select
+        end associate
+
+        if (present(interval)) then
+            associate ( a => interval(1) , b => interval(2) )
+            x = 0.5*(b-a)*x+0.5*(b+a)
+            w = 0.5*(b-a)*w
+            end associate
+        end if
+    end subroutine
+
+    pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
+        real(dp), intent(out) :: x(:), w(:)
+        real(dp), intent(in), optional :: interval(2)
+
+        associate (N => size(x)-1)
+        select case (N)
+            case (1)
+                x = [-1._dp, 1._dp]
+                w = [ 1._dp, 1._dp]
+            case default
+                block
+                integer :: i,j
+                real(dp) :: leg, dleg, delta
+
+                x(1)   = -1._dp
+                x(N+1) =  1._dp
+                w(1)   =  2._dp/(N*(N+1._dp))
+                w(N+1) =  2._dp/(N*(N+1._dp))
+
+                do i = 1, int(floor((N+1)/2._dp)-1)
+                    x(i+1) = -cos( (i+0.25_dp)*PI/N  - 3/(8*N*PI*(i+0.25_dp)))
+                    do j = 0, newton_iters-1
+                        leg  = legendre(N+1,x(i+1)) - legendre(N-1,x(i+1))
+                        dleg = dlegendre(N+1,x(i+1)) - dlegendre(N-1,x(i+1))
+                        delta = -leg/dleg
+                        x(i+1) = x(i+1) + delta
+                        if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
+                    end do
+                    x(N-i+1) = -x(i+1)
+
+                    leg = legendre(N, x(i+1))
+                    w(i+1)   = 2._dp/(N*(N+1._dp)*leg**2) 
+                    w(N-i+1) = w(i+1)
+                end do
+
+                if (mod(N,2) == 0) then
+                    x(N/2+1) = 0.0
+
+                    leg = legendre(N, 0.0_dp)
+                    w(N/2+1)   = 2._dp/(N*(N+1._dp)*leg**2) 
+                end if
+                end block
+        end select
+        end associate
+        
+        if (present(interval)) then
+            associate ( a => interval(1) , b => interval(2) )
+            x = 0.5*(b-a)*x+0.5*(b+a)
+            w = 0.5*(b-a)*w
+            end associate
+        end if
+    end subroutine
+end submodule