diff --git a/Modules/ThirdParty/VNL/src/CMakeLists.txt b/Modules/ThirdParty/VNL/src/CMakeLists.txt
index 7a7cb9af078..c51bc55b91c 100644
--- a/Modules/ThirdParty/VNL/src/CMakeLists.txt
+++ b/Modules/ThirdParty/VNL/src/CMakeLists.txt
@@ -26,7 +26,6 @@ endforeach()
 
 foreach(
   exe
-  netlib_integral_test
   netlib_lbfgs_example
   netlib_lbfgsb_example
   netlib_slamch_test
diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_integral.cxx b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_integral.cxx
index 1eec4eaf05a..ad421e55c63 100644
--- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_integral.cxx
+++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_integral.cxx
@@ -1,3 +1,4 @@
+#define ITK_LEGACY_TEST
 // not used? #include <iostream>
 #include <cmath>
 #include <utility>
diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_adaptsimpson_integral.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_adaptsimpson_integral.h
index 0bec1f448c6..a37a69e8184 100644
--- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_adaptsimpson_integral.h
+++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_adaptsimpson_integral.h
@@ -1,5 +1,18 @@
 #ifndef VNL_ADAPTSIMPSON_INTEGRAL_H_
 #define VNL_ADAPTSIMPSON_INTEGRAL_H_
+
+#if __has_include(<itkConfigure.h>)
+#  include <itkConfigure.h>
+#  if defined(ITK_FUTURE_LEGACY_REMOVE)
+#    error "vnl_adaptsimpson_integral was removed; supply an adaptive Simpson rule directly."
+#  elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST)
+#    if defined(_MSC_VER)
+#      pragma message("vnl_adaptsimpson_integral is deprecated; supply an adaptive Simpson rule directly.")
+#    else
+#      warning "vnl_adaptsimpson_integral is deprecated; supply an adaptive Simpson rule directly."
+#    endif
+#  endif
+#endif
 //:
 // \file
 // \author Kongbin Kang at Brown
diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.cxx b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.cxx
index 722a07840e7..0ed1db6fa1c 100644
--- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.cxx
+++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.cxx
@@ -1,5 +1,4 @@
 #include "vnl_simpson_integral.h"
-#include <vnl/algo/vnl_netlib.h>
 
 double
 vnl_simpson_integral::int_fnct_(double * x)
@@ -10,13 +9,25 @@ vnl_simpson_integral::int_fnct_(double * x)
 double
 vnl_simpson_integral::integral(vnl_integrant_fnct * f, double a, double b, long n)
 {
-
-  double res = 0;
-
   // set the function
   pfnct_ = f;
 
-  v3p_netlib_simpru_(&vnl_simpson_integral::int_fnct_, &a, &b, &n, &res);
+  // Composite Simpson rule (Mathews Algorithm 7.2): n intervals, 2n subintervals.
+  const double h = (b - a) / (2 * n);
+
+  double sumeven = 0.0;
+  for (long k = 1; k <= n - 1; ++k)
+  {
+    double x = a + h * 2 * k;
+    sumeven += int_fnct_(&x);
+  }
+
+  double sumodd = 0.0;
+  for (long k = 1; k <= n; ++k)
+  {
+    double x = a + h * (2 * k - 1);
+    sumodd += int_fnct_(&x);
+  }
 
-  return res;
+  return h * (int_fnct_(&a) + int_fnct_(&b) + 2 * sumeven + 4 * sumodd) / 3;
 }
diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.h
index 8af81e1aab7..5018a94c034 100644
--- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.h
+++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.h
@@ -1,5 +1,18 @@
 #ifndef VNL_SIMPSON_INTEGRAL_H_
 #define VNL_SIMPSON_INTEGRAL_H_
+
+#if __has_include(<itkConfigure.h>)
+#  include <itkConfigure.h>
+#  if defined(ITK_FUTURE_LEGACY_REMOVE)
+#    error "vnl_simpson_integral was removed; supply a composite Simpson rule (Mathews Algorithm 7.2) directly."
+#  elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST)
+#    if defined(_MSC_VER)
+#      pragma message("vnl_simpson_integral is deprecated; supply a composite Simpson rule directly.")
+#    else
+#      warning "vnl_simpson_integral is deprecated; supply a composite Simpson rule directly."
+#    endif
+#  endif
+#endif
 //:
 // \file
 // \author Kongbin Kang at Brown
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt
index 1c8ced056cf..65ef65ea79d 100644
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt
+++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt
@@ -166,10 +166,6 @@ set(V3P_NETLIB_linalg_SOURCES
      linalg/lsmrBase.cxx linalg/lsmrBase.h
      linalg/lsmrDense.cxx linalg/lsmrDense.h
    )
-set(V3P_NETLIB_mathews_SOURCES
-     mathews/simpson.c mathews/simpson.h
-     mathews/trapezod.c mathews/trapezod.h
-   )
 set(V3P_NETLIB_sparse_SOURCES
      sparse/spAllocate.c
      sparse/spBuild.c
@@ -192,7 +188,6 @@ set(v3p_netlib_sources
      ${V3P_NETLIB_napack_SOURCES}
      ${V3P_NETLIB_opt_SOURCES}
      ${V3P_NETLIB_linalg_SOURCES}
-     ${V3P_NETLIB_mathews_SOURCES}
      ${V3P_NETLIB_sparse_SOURCES}
    )
 
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/README.md b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/README.md
index 4c711428a7a..6afa808eca4 100644
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/README.md
+++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/README.md
@@ -73,7 +73,7 @@ respecitve sources.
 5. Convert the sources to C using `f2c` and replace the `f2c.h` header
 inclusion with `v3p_netlib.h`:
     ```bash
-    for d in blas linpack temperton eispack laso arpack lapack/complex16 lapack/double lapack/single lapack/util napack minpack opt linalg toms datapac mathews; do
+    for d in blas linpack temperton eispack laso arpack lapack/complex16 lapack/double lapack/single lapack/util napack minpack opt linalg toms datapac; do
       for f in ${d}/*.f; do
         b=`echo "$f" | sed 's/.f$//'`
         if [ ! -f "${b}.c" ]; then
@@ -120,7 +120,7 @@ be left alone.
 with the mangled interface.  Include them all in `v3p_netlib_prototypes.h`.
     ```bash
     echo "/* Include prototype headers. */" > v3p_netlib_prototypes.h
-    for f in blas/*.P linpack/*.P temperton/*.P eispack/*.P laso/*.P arpack/*.P lapack/*/*.P napack/*.P minpack/*.P opt/*.P linalg/*.P toms/*.P datapac/*.P mathews/*.P; do
+    for f in blas/*.P linpack/*.P temperton/*.P eispack/*.P laso/*.P arpack/*.P lapack/*/*.P napack/*.P minpack/*.P opt/*.P linalg/*.P toms/*.P datapac/*.P; do
       b=`echo "$f" | sed 's/.P//'`
       if [ ! -f "${b}.h" ] ; then
         echo "Converting prototype $b"
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/adaquad.f b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/adaquad.f
deleted file mode 100644
index 3477cf0b0fc..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/adaquad.f
+++ /dev/null
@@ -1,131 +0,0 @@
-C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
-C     To accompany the text:
-C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
-C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
-C     This free software is complements of the author.
-C
-C     Algorithm 7.5 (Adaptive Quadrature Using Simpson's Rule).
-C     Section 7.4, Adaptive Quadrature, Page 389
-C
-
-C     add missing variable F in Refine subrutine.
-      SUBROUTINE AdaptQuad(F,A,B,Tol,SRmat,Integral,ErrBdd,M,State)
-      INTEGER M,State
-      DOUBLE PRECISION A,B,Tol,SRmat,Integral,ErrBdd
-      INTEGER J,K,N,Iterating,Done
-      DOUBLE PRECISION Sum1,Sum2,SRvec
-      DIMENSION SRmat(1:101,1:11),SRvec(1:11)
-      EXTERNAL F
-      Iterating = 0
-      Done = 1
-      CALL Srule(F,A,B,Tol,SRvec)
-      DO K=1,11
-        SRmat(1, K) = SRvec(K)
-      ENDDO
-      M = 1
-      State = Iterating
-      DO WHILE (State .EQ. Iterating)
-        N = M
-        DO J=N,1,-1
-          CALL Refine(F,J,SRmat,M,State)
-        ENDDO
-      ENDDO
-      Sum1 = 0
-      Sum2 = 0
-      DO J=1,M
-        Sum1 = Sum1 + SRmat(J, 8)
-        Sum2 = Sum2 + Abs(SRmat(J, 9))
-      ENDDO
-      Integral = Sum1
-      ErrBdd = Sum2
-      RETURN
-      END
-
-      SUBROUTINE Refine(F, P,SRmat,M,State)
-      INTEGER P,M,State
-      DOUBLE PRECISION SRmat
-      INTEGER J,K,Iterating,Done
-      DOUBLE PRECISION A,B,C,Err,Fa,Fb,Fc,S,S2,Tol,Tol2,Err,Check
-      DOUBLE PRECISION SR0vec,SR1vec,SR2vec
-      DIMENSION SRmat(1:101,1:11)
-      DIMENSION SR0vec(1:11),SR1vec(1:11),SR2vec(1:11)
-      EXTERNAL F
-      Iterating = 0
-      Done = 1
-      State = Done
-      DO K=1,11
-        SR0vec(K) = SRmat(P, K)
-      ENDDO
-      A = SR0vec(1)
-      C = SR0vec(2)
-      B = SR0vec(3)
-      Fa = SR0vec(4)
-      Fc = SR0vec(5)
-      Fb = SR0vec(6)
-      S = SR0vec(7)
-      S2 = SR0vec(8)
-      Err = SR0vec(9)
-      Tol = SR0vec(10)
-      Check = SR0vec(11)
-      IF (Check .EQ. 1) RETURN
-      Tol2 = Tol / 2
-      CALL Srule(F, A, C, Tol2, SR1vec)
-      CALL Srule(F, C, B, Tol2, SR2vec)
-      Err = ABS(SR0vec(7) - SR1vec(7) - SR2vec(7)) / 10
-      IF (Err .LT. Tol) THEN
-        SR0vec(11) = 1
-      ENDIF
-      IF (Err .LT. Tol) THEN
-        DO K=1,11
-          SRmat(P, K) = SR0vec(K)
-        ENDDO
-        SRmat(P, 8) = SR1vec(7) + SR2vec(7)
-        SRmat(P, 9) = Err
-      ELSE
-        DO J=(M + 1),P,-1
-          DO K=1,11
-            SRmat(J, K) = SRmat(J - 1, K)
-          ENDDO
-        ENDDO
-        M = M + 1
-        DO K=1,11
-          SRmat(P, K) = SR1vec(K)
-        ENDDO
-        DO K=1,11
-          SRmat(P + 1, K) = SR2vec(K)
-        ENDDO
-        State = Iterating
-      ENDIF
-      RETURN
-      END
-
-      SUBROUTINE Srule(F,A,B,Tol0,SRvec)
-      DOUBLE PRECISION A,B,Tol0,SRvec
-      DOUBLE PRECISION C,H,Fa,Fb,Fc,S,S2,Tol1,Err,Check
-      DIMENSION SRvec(1:11)
-      EXTERNAL F
-        H = (B - A) / 2
-        C = (A + B) / 2
-        Fa = F(A)
-        Fc = F(C)
-        Fb = F(B)
-        S = H * (F(A) + 4 * F(C) + F(B)) / 3
-        S2 = S
-        Tol1 = Tol0
-        Err = Tol0
-        Check = 0
-        SRvec(1) = A
-        SRvec(2) = C
-        SRvec(3) = B
-        SRvec(4) = Fa
-        SRvec(5) = Fc
-        SRvec(6) = Fb
-        SRvec(7) = S
-        SRvec(8) = S2
-        SRvec(9) = Err
-        SRvec(10) = Tol1
-        SRvec(11) = Check
-      RETURN
-      END
-
-
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.c b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.c
deleted file mode 100644
index b72f6932c2d..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.c
+++ /dev/null
@@ -1,136 +0,0 @@
-/* mathews/simpson.f -- translated by f2c (version 20050501).
-   You must link the resulting object file with libf2c:
-        on Microsoft Windows system, link with libf2c.lib;
-        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
-        or, if you install libf2c.a in a standard place, with -lf2c -lm
-        -- in that order, at the end of the command line, as in
-                cc *.o -lf2c -lm
-        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
-
-                http://www.netlib.org/f2c/libf2c.zip
-*/
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-#include "v3p_netlib.h"
-
-/*     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 */
-/*     To accompany the text: */
-/*     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 */
-/*     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. */
-/*     This free software is complements of the author. */
-
-/*     Algorithm 7.2 (Composite Simpson Rule). */
-/*     Section 7.2, Composite Trapezoidal and Simpson's Rule, Page 365 */
-
-/*     comment added by Kongbin Kang */
-/*     F: integrand function */
-/*     A: lower integration limit */
-/*     B: higher integration limit */
-/*     M: number of intervals. Notice, the subintervals used is 2M */
-/*     Srule: output parameter to store simpson rule result */
-/*<       SUBROUTINE SIMPRU(F,A,B,M,Srule) >*/
-/* Subroutine */ int simpru_(doublereal (*f)(doublereal*),
-                             doublereal *a, doublereal *b, integer *m,
-                             doublereal *srule)
-{
-    /* System generated locals */
-    integer i__1;
-
-    /* Local variables */
-    doublereal h__;
-    integer k;
-    doublereal x, sum, sumodd, sumeven;
-
-/*<       INTEGER K,M >*/
-/*<       DOUBLE PRECISION A,B,H,Sum,SumEven,SumOdd,Srule,X >*/
-/*<       EXTERNAL F >*/
-/*<       H=(B-A)/(2*M) >*/
-    h__ = (*b - *a) / (*m << 1);
-/*<       SumEven=0 >*/
-    sumeven = 0.;
-/*<       DO K=1,(M-1) >*/
-    i__1 = *m - 1;
-    for (k = 1; k <= i__1; ++k) {
-/*<         X=A+H*2*K >*/
-        x = *a + h__ * 2 * k;
-/*<         SumEven=SumEven+F(X) >*/
-        sumeven += (*f)(&x);
-/*<       ENDDO >*/
-    }
-/*<       SumOdd=0 >*/
-    sumodd = 0.;
-/*<       DO K=1,M >*/
-    i__1 = *m;
-    for (k = 1; k <= i__1; ++k) {
-/*<         X=A+H*(2*K-1) >*/
-        x = *a + h__ * ((k << 1) - 1);
-/*<         SumOdd=SumOdd+F(X) >*/
-        sumodd += (*f)(&x);
-/*<       ENDDO >*/
-    }
-/*<       Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3 >*/
-    sum = h__ * ((*f)(a) + (*f)(b) + sumeven * 2 + sumodd * 4) / 3;
-/*<       Srule=Sum >*/
-    *srule = sum;
-/*<       RETURN >*/
-    return 0;
-/*<       END >*/
-} /* simpru_ */
-
-/*<       SUBROUTINE XSIMPRU(F,A,B,M,Srule) >*/
-/* Subroutine */ int xsimpru_(doublereal (*f)(doublereal*),
-                              doublereal *a, doublereal *b, integer *
-                              m, doublereal *srule)
-{
-    /* System generated locals */
-    integer i__1;
-
-    /* Local variables */
-    doublereal h__;
-    integer k;
-    doublereal x, sum, sumodd, sumeven;
-
-/*     This subroutine uses labeled DO loop(s). */
-/*<       INTEGER K,M >*/
-/*<       DOUBLE PRECISION A,B,H,Sum,SumEven,SumOdd,Srule,X >*/
-/*<       EXTERNAL F >*/
-/*<       H=(B-A)/(2*M) >*/
-    h__ = (*b - *a) / (*m << 1);
-/*<       SumEven=0 >*/
-    sumeven = 0.;
-/*<       DO 10 K=1,(M-1) >*/
-    i__1 = *m - 1;
-    for (k = 1; k <= i__1; ++k) {
-/*<         X=A+H*2*K >*/
-        x = *a + h__ * 2 * k;
-/*<         SumEven=SumEven+F(X) >*/
-        sumeven += (*f)(&x);
-/*< 10    CONTINUE >*/
-/* L10: */
-    }
-/*<       SumOdd=0 >*/
-    sumodd = 0.;
-/*<       DO 20 K=1,M >*/
-    i__1 = *m;
-    for (k = 1; k <= i__1; ++k) {
-/*<         X=A+H*(2*K-1) >*/
-        x = *a + h__ * ((k << 1) - 1);
-/*<         SumOdd=SumOdd+F(X) >*/
-        sumodd += (*f)(&x);
-/*< 20    CONTINUE >*/
-/* L20: */
-    }
-/*<       Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3 >*/
-    sum = h__ * ((*f)(a) + (*f)(b) + sumeven * 2 + sumodd * 4) / 3;
-/*<       Srule=Sum >*/
-    *srule = sum;
-/*<       RETURN >*/
-    return 0;
-/*<       END >*/
-} /* xsimpru_ */
-
-#ifdef __cplusplus
-        }
-#endif
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.f b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.f
deleted file mode 100644
index 59063f55fb6..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.f
+++ /dev/null
@@ -1,59 +0,0 @@
-C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
-C     To accompany the text:
-C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
-C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
-C     This free software is complements of the author.
-C
-C     Algorithm 7.2 (Composite Simpson Rule).
-C     Section 7.2, Composite Trapezoidal and Simpson's Rule, Page 365
-C
-
-C     comment added by Kongbin Kang
-C     F: integrand function
-C     A: lower integration limit
-C     B: higher integration limit
-C     M: number of intervals. Notice, the subintervals used is 2M
-C     Srule: output parameter to store simpson rule result
-
-      SUBROUTINE SIMPRU(F,A,B,M,Srule)
-      INTEGER K,M
-      DOUBLE PRECISION A,B,H,Sum,SumEven,SumOdd,Srule,X
-      EXTERNAL F
-      H=(B-A)/(2*M)
-      SumEven=0
-      DO K=1,(M-1)
-        X=A+H*2*K
-        SumEven=SumEven+F(X)
-      ENDDO
-      SumOdd=0
-      DO K=1,M
-        X=A+H*(2*K-1)
-        SumOdd=SumOdd+F(X)
-      ENDDO
-      Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3
-      Srule=Sum
-      RETURN
-      END
-
-      SUBROUTINE XSIMPRU(F,A,B,M,Srule)
-C     This subroutine uses labeled DO loop(s).
-      INTEGER K,M
-      DOUBLE PRECISION A,B,H,Sum,SumEven,SumOdd,Srule,X
-      EXTERNAL F
-      H=(B-A)/(2*M)
-      SumEven=0
-      DO 10 K=1,(M-1)
-        X=A+H*2*K
-        SumEven=SumEven+F(X)
-10    CONTINUE
-      SumOdd=0
-      DO 20 K=1,M
-        X=A+H*(2*K-1)
-        SumOdd=SumOdd+F(X)
-20    CONTINUE
-      Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3
-      Srule=Sum
-      RETURN
-      END
-
-
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.h b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.h
deleted file mode 100644
index e5ec0e8546c..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/simpson.h
+++ /dev/null
@@ -1,14 +0,0 @@
-extern int v3p_netlib_simpru_(
-  v3p_netlib_doublereal (*f)(v3p_netlib_doublereal*),
-  v3p_netlib_doublereal *a,
-  v3p_netlib_doublereal *b,
-  v3p_netlib_integer *m,
-  v3p_netlib_doublereal *srule
-  );
-extern int v3p_netlib_xsimpru_(
-  v3p_netlib_doublereal (*f)(v3p_netlib_doublereal*),
-  v3p_netlib_doublereal *a,
-  v3p_netlib_doublereal *b,
-  v3p_netlib_integer *m,
-  v3p_netlib_doublereal *srule
-  );
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.c b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.c
deleted file mode 100644
index 87d490add5b..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.c
+++ /dev/null
@@ -1,109 +0,0 @@
-/* mathews/trapezod.f -- translated by f2c (version 20050501).
-   You must link the resulting object file with libf2c:
-        on Microsoft Windows system, link with libf2c.lib;
-        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
-        or, if you install libf2c.a in a standard place, with -lf2c -lm
-        -- in that order, at the end of the command line, as in
-                cc *.o -lf2c -lm
-        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
-
-                http://www.netlib.org/f2c/libf2c.zip
-*/
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-#include "v3p_netlib.h"
-
-/*     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 */
-/*     To accompany the text: */
-/*     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 */
-/*     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. */
-/*     This free software is complements of the author. */
-
-/*     Algorithm 7.1 (Composite Trapezoidal Rule). */
-/*     Section 7.2, Composite Trapezoidal and Simpson's Rule, Page 365 */
-
-/*<       SUBROUTINE TRAPRU(F,A,B,M,Trule) >*/
-/* Subroutine */ int trapru_(
-  v3p_netlib_doublereal (*f)(v3p_netlib_doublereal*),
-  doublereal *a, doublereal *b, integer *m, doublereal *trule
-  )
-{
-    /* System generated locals */
-    integer i__1;
-
-    /* Local variables */
-    doublereal h__;
-    integer k;
-    doublereal x, sum;
-
-/*<       INTEGER K,M >*/
-/*<       DOUBLE PRECISION A,B,H,Sum,Trule,X >*/
-/*<       EXTERNAL F >*/
-/*<       H=(B-A)/M >*/
-    h__ = (*b - *a) / *m;
-/*<       Sum=0 >*/
-    sum = 0.;
-/*<       DO K=1,M-1 >*/
-    i__1 = *m - 1;
-    for (k = 1; k <= i__1; ++k) {
-/*<         X=A+H*K >*/
-        x = *a + h__ * k;
-/*<         Sum=Sum+F(X) >*/
-        sum += (*f)(&x);
-/*<       ENDDO >*/
-    }
-/*<       Sum=H*(F(A)+F(B)+2*Sum)/2 >*/
-    sum = h__ * ((*f)(a) + (*f)(b) + sum * 2) / 2;
-/*<       Trule=Sum >*/
-    *trule = sum;
-/*<       RETURN >*/
-    return 0;
-/*<       END >*/
-} /* trapru_ */
-
-/*<       SUBROUTINE XTRAPRU(F,A,B,M,Trule) >*/
-/* Subroutine */ int xtrapru_(
-  v3p_netlib_doublereal (*f)(v3p_netlib_doublereal*),
-  doublereal *a, doublereal *b, integer * m, doublereal *trule
-  )
-{
-    /* System generated locals */
-    integer i__1;
-
-    /* Local variables */
-    doublereal h__;
-    integer k;
-    doublereal x, sum;
-
-/*     This subroutine uses labeled DO loop(s). */
-/*<       INTEGER K,M >*/
-/*<       DOUBLE PRECISION A,B,H,Sum,Trule,X >*/
-/*<       EXTERNAL F >*/
-/*<       H=(B-A)/M >*/
-    h__ = (*b - *a) / *m;
-/*<       Sum=0 >*/
-    sum = 0.;
-/*<       DO 10 K=1,M-1 >*/
-    i__1 = *m - 1;
-    for (k = 1; k <= i__1; ++k) {
-/*<         X=A+H*K >*/
-        x = *a + h__ * k;
-/*<         Sum=Sum+F(X) >*/
-        sum += (*f)(&x);
-/*< 10    CONTINUE >*/
-/* L10: */
-    }
-/*<       Sum=H*(F(A)+F(B)+2*Sum)/2 >*/
-    sum = h__ * ((*f)(a) + (*f)(b) + sum * 2) / 2;
-/*<       Trule=Sum >*/
-    *trule = sum;
-/*<       RETURN >*/
-    return 0;
-/*<       END >*/
-} /* xtrapru_ */
-
-#ifdef __cplusplus
-        }
-#endif
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.f b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.f
deleted file mode 100644
index 8122daff934..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.f
+++ /dev/null
@@ -1,47 +0,0 @@
-C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
-C     To accompany the text:
-C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
-C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
-C     This free software is complements of the author.
-C
-C     Algorithm 7.1 (Composite Trapezoidal Rule).
-C     Section 7.2, Composite Trapezoidal and Simpson's Rule, Page 365
-C
-      SUBROUTINE TRAPRU(F,A,B,M,Trule)
-      INTEGER K,M
-      DOUBLE PRECISION A,B,H,Sum,Trule,X
-      EXTERNAL F
-      H=(B-A)/M
-      Sum=0
-      DO K=1,M-1
-        X=A+H*K
-        Sum=Sum+F(X)
-      ENDDO
-      Sum=H*(F(A)+F(B)+2*Sum)/2
-      Trule=Sum
-      RETURN
-      END
-
-      SUBROUTINE XTRAPRU(F,A,B,M,Trule)
-C     This subroutine uses labeled DO loop(s).
-      INTEGER K,M
-      DOUBLE PRECISION A,B,H,Sum,Trule,X
-      EXTERNAL F
-      H=(B-A)/M
-      Sum=0
-      DO 10 K=1,M-1
-        X=A+H*K
-        Sum=Sum+F(X)
-10    CONTINUE
-      Sum=H*(F(A)+F(B)+2*Sum)/2
-      Trule=Sum
-      RETURN
-      END
-
-
-
-
-
-
-
-
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.h b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.h
deleted file mode 100644
index 8af67b09583..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/mathews/trapezod.h
+++ /dev/null
@@ -1,14 +0,0 @@
-extern int v3p_netlib_trapru_(
-  v3p_netlib_doublereal (*f)(v3p_netlib_doublereal*),
-  v3p_netlib_doublereal *a,
-  v3p_netlib_doublereal *b,
-  v3p_netlib_integer *m,
-  v3p_netlib_doublereal *trule
-  );
-extern int v3p_netlib_xtrapru_(
-  v3p_netlib_doublereal (*f)(v3p_netlib_doublereal*),
-  v3p_netlib_doublereal *a,
-  v3p_netlib_doublereal *b,
-  v3p_netlib_integer *m,
-  v3p_netlib_doublereal *trule
-  );
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/CMakeLists.txt b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/CMakeLists.txt
index 464e277b7d0..3d70fa29532 100644
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/CMakeLists.txt
+++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/CMakeLists.txt
@@ -23,8 +23,4 @@ if(BUILD_TESTING)
   add_executable( netlib_slamch_test slamch-test.c )
   add_test( NAME netlib_test_slamch COMMAND netlib_slamch_test )
   target_link_libraries(netlib_slamch_test itkv3p_netlib)
-  # test
-  add_executable( netlib_integral_test integral-test.c )
-  add_test( NAME netlib_test_integral COMMAND netlib_integral_test )
-  target_link_libraries(netlib_integral_test itkv3p_netlib)
 endif()
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/integral-test.c b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/integral-test.c
deleted file mode 100644
index 17b778ce95c..00000000000
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/tests/integral-test.c
+++ /dev/null
@@ -1,39 +0,0 @@
-#include "v3p_netlib.h"
-
-#include <stdio.h>
-
-double f(double *x)
-{
-  return (*x)/(1+(*x)*(*x));
-}
-
-void test_simpson_integral()
-{
-  double a = 0;
-  double b = 1;
-  double res;
-  long n = 100;
-
-  v3p_netlib_simpru_(&f, &a, &b, &n, &res);
-  printf("simpson integral of x/(1+x^2) from 0 to 1 (%ld grids) is %2.10f\n", n, res);
-}
-
-void test_trapezod_integral()
-{
-  double a = 0;
-  double b = 1;
-  double res;
-  long n = 500;
-
-  v3p_netlib_trapru_(&f, &a, &b, &n, &res);
-  printf("trapezod integral of x/(1+x^2) from 0 to 1 (%ld grids) is %f\n", n, res);
-}
-
-
-
-int main()
-{
-  test_simpson_integral();
-  test_trapezod_integral();
-  return 0;
-}
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_mangle.h b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_mangle.h
index dcb30b83363..c8484ad3f29 100644
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_mangle.h
+++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_mangle.h
@@ -169,7 +169,6 @@ nm libv3p_netlib.a |grep " [TR] "
 #define sgbmv_ v3p_netlib_sgbmv_
 #define sgemv_ v3p_netlib_sgemv_
 #define sger_ v3p_netlib_sger_
-#define simpru_ v3p_netlib_simpru_
 #define slamc1_ v3p_netlib_slamc1_
 #define slamc2_ v3p_netlib_slamc2_
 #define slamc3_ v3p_netlib_slamc3_
@@ -206,13 +205,10 @@ nm libv3p_netlib.a |grep " [TR] "
 #define tql1_ v3p_netlib_tql1_
 #define tql2_ v3p_netlib_tql2_
 #define tqlrat_ v3p_netlib_tqlrat_
-#define trapru_ v3p_netlib_trapru_
 #define tred1_ v3p_netlib_tred1_
 #define tred2_ v3p_netlib_tred2_
 #define urand_ v3p_netlib_urand_
 #define xerbla_ v3p_netlib_xerbla_
-#define xsimpru_ v3p_netlib_xsimpru_
-#define xtrapru_ v3p_netlib_xtrapru_
 #define zaxpy_ v3p_netlib_zaxpy_
 #define zcopy_ v3p_netlib_zcopy_
 #define zdotc_ v3p_netlib_zdotc_
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_prototypes.h b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_prototypes.h
index 0177f3de64b..e2bed118876 100644
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_prototypes.h
+++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_prototypes.h
@@ -94,5 +94,3 @@
 #include "napack/cg.h"
 #include "opt/lbfgs.h"
 #include "opt/lbfgsb.h"
-#include "mathews/simpson.h"
-#include "mathews/trapezod.h"
diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_unmangle.h b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_unmangle.h
index ce1f2463d40..3145e45b8ba 100644
--- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_unmangle.h
+++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/v3p_netlib_unmangle.h
@@ -160,7 +160,6 @@
 #undef sgbmv_
 #undef sgemv_
 #undef sger_
-#undef simpru_
 #undef slamc1_
 #undef slamc2_
 #undef slamc3_
@@ -197,13 +196,10 @@
 #undef tql1_
 #undef tql2_
 #undef tqlrat_
-#undef trapru_
 #undef tred1_
 #undef tred2_
 #undef urand_
 #undef xerbla_
-#undef xsimpru_
-#undef xtrapru_
 #undef zaxpy_
 #undef zcopy_
 #undef zdotc_