Merge pull request #660 from hsnyder/master

Fix erroneous gaussian quadrature points in gauss_legendre

Author: milancurcic

src/stdlib_quadrature_gauss.f90

@@ -56,8 +56,6 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
         if (present(interval)) then
             associate ( a => interval(1) , b => interval(2) )
             x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
-            x(1)       = interval(1)
-            x(size(x)) = interval(2)
             w = 0.5_dp*(b-a)*w
             end associate
         end if

src/tests/quadrature/test_gauss.f90

@@ -21,7 +21,8 @@ subroutine collect_gauss(testsuite)
             new_unittest("gauss-lobatto-analytic", test_gauss_lobatto_analytic), &
             new_unittest("gauss-lobatto-5", test_gauss_lobatto_5), &
             new_unittest("gauss-lobatto-32", test_gauss_lobatto_32), &
-            new_unittest("gauss-lobatto-64", test_gauss_lobatto_64) &
+            new_unittest("gauss-lobatto-64", test_gauss_lobatto_64), &
+            new_unittest("gauss-github-issue-619", test_fix_github_issue619) &
             ]
     end subroutine
 
@@ -48,6 +49,25 @@ subroutine test_gauss_analytic(error)
 
     end subroutine
 
+    subroutine test_fix_github_issue619(error)
+        !> See github issue https://github.com/fortran-lang/stdlib/issues/619
+        type(error_type), allocatable, intent(out) :: error
+        integer :: i
+
+        ! test the values of nodes and weights
+        i = 5
+        block 
+            real(dp), dimension(i) :: x1,w1,x2,w2
+            call gauss_legendre(x1,w1)
+            call gauss_legendre(x2,w2,interval=[-1._dp, 1._dp])
+
+            call check(error, all(abs(x1-x2) < 2*epsilon(x1(1))))
+            if (allocated(error)) return
+            call check(error, all(abs(w1-w2) < 2*epsilon(w1(1))))
+        end block
+
+    end subroutine
+
     subroutine test_gauss_5(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error