Improved GaussLegendre and LegendreFactory

Now, the computation of the nodes and weights of the quadrature rules are done using FastGL (https://sourceforge.net/projects/fastgausslegendrequadrature/), one of the fastest and most accurate algorithm to do so.

CMakeLists.txt

@@ -677,6 +677,7 @@ install (FILES COPYING
                 COPYING.ev3
                 COPYING.exprtk
                 COPYING.faddeeva
+                COPYING.fastgl
                 COPYING.kendall
                 COPYING.kissfft
                 COPYING.cephes

COPYING.fastgl

@@ -0,0 +1,20 @@
+//*******************************************
+//   Copyright (C) 2014 by Ignace Bogaert   *
+//*******************************************
+
+// This software package is based on the paper
+//    I. Bogaert, "Iteration-Free Computation of Gauss-Legendre Quadrature Nodes and Weights",
+//    to be published in the SIAM Journal of Scientific Computing.
+
+// The main features of this software are:
+// - Speed: due to the simple formulas and the O(1) complexity computation of individual Gauss-Legendre 
+//   quadrature nodes and weights. This makes it compatible with parallel computing paradigms.
+// - Accuracy: the error on the nodes and weights is within a few ulps (see the paper for details).
+
+// Disclaimer:
+// THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR 
+// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
+// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

ChangeLog

@@ -15,6 +15,8 @@
 === Python module ===
 
 === Miscellaneous ===
+ * New Gauss Legendre algorithm to generate weights and nodes, relying on
+	FastGL (https://sourceforge.net/projects/fastgausslegendrequadrature/)
 
 
 

LICENSE

@@ -10,3 +10,4 @@ This library bundles several third-party codes with various licenses compatible
 - KissFFT (lib/src/Base/Algo/kissfft.hh), under BSD license, see COPYING.kissfft
 - KS distribution from Cephes library (lib/src/Uncertainty/Distribution/cephes/*) under BSD license, see COPYING.cephes
 - TNC optimization solver (lib/src/Base/Optim/algotnc.*) under Expat license, see COPYING.tnc
+- Gauss Legendre quadrature from FastGL library (lib/src/Base/Algo/fastgl*) under BSD license, see COPYING.fastGL

lib/src/Base/Algo/CMakeLists.txt

@@ -36,6 +36,7 @@ ot_add_source_file (IntegrationAlgorithm.cxx)
 ot_add_source_file (FilonQuadrature.cxx)
 ot_add_source_file (GaussKronrod.cxx)
 ot_add_source_file (GaussKronrodRule.cxx)
+ot_add_source_file (fastgl.cpp)
 ot_add_source_file (GaussLegendre.cxx)
 ot_add_source_file (FejerAlgorithm.cxx)
 ot_add_source_file (IteratedQuadrature.cxx)

lib/src/Base/Algo/GaussLegendre.cxx

@@ -22,7 +22,7 @@
 #include "openturns/Tuples.hxx"
 #include "openturns/Exception.hxx"
 #include "openturns/PersistentObjectFactory.hxx"
-#include "openturns/Lapack.hxx"
+#include "fastgl.h"
 
 BEGIN_NAMESPACE_OPENTURNS
 
@@ -106,45 +106,14 @@ void GaussLegendre::generateNodesAndWeights()
     if (!(integrationNodesNumber > 0)) throw InvalidArgumentException(HERE) << "Error: the discretization must be positive, here discretization[" << i << "] has " << integrationNodesNumber << "nodes.";
     // Check if we already computed this 1D rule
     // We use the value 'dimension' as a guard
-    UnsignedInteger indexAlreadyComputed = dimension;
-    for (UnsignedInteger j = 0; j < i; ++j)
-      if (discretization_[j] == integrationNodesNumber)
+    marginalNodes[i] = Point(integrationNodesNumber);
+    marginalWeights[i] = Point(integrationNodesNumber);
+    for (UnsignedInteger j = 0; j < integrationNodesNumber; ++j)
       {
-        indexAlreadyComputed = j;
-        break;
-      }
-    // If indexAlreadyComputed < dimension we found a match
-    if (indexAlreadyComputed < dimension)
-    {
-      marginalNodes[i] = marginalNodes[indexAlreadyComputed];
-      marginalWeights[i] = marginalWeights[indexAlreadyComputed];
-    } // A match found
-    else
-    {
-      marginalNodes[i] = Point(integrationNodesNumber);
-      marginalWeights[i] = Point(integrationNodesNumber);
-      // First, build a symmetric tridiagonal matrix whose eigenvalues are the nodes of the
-      // gauss integration rule
-      char jobz('V');
-      int ljobz(1);
-      Point d(integrationNodesNumber);
-      Point e(integrationNodesNumber);
-      for (UnsignedInteger k = 1; k < integrationNodesNumber; ++k) e[k - 1] = 0.5 / std::sqrt(1.0 - std::pow(2.0 * k, -2));
-      int ldz(integrationNodesNumber);
-      SquareMatrix z(integrationNodesNumber);
-      Point work(2 * integrationNodesNumber - 2);
-      int info;
-      int n = static_cast<int>(integrationNodesNumber);
-      dstev_(&jobz, &n, &d[0], &e[0], &z(0, 0), &ldz, &work[0], &info, &ljobz);
-      if (info != 0) throw InternalException(HERE) << "Lapack DSTEV: error code=" << info;
-      for (UnsignedInteger k = 0; k < integrationNodesNumber; ++k)
-      {
-        // Nodes
-        marginalNodes[i][k] = 0.5 * (1.0 + d[k]);
-        // Weights
-        marginalWeights[i][k] = std::pow(z(0, k), 2);
-      }
-    } // No match found
+	fastgl::QuadPair p(fastgl::GLPair(integrationNodesNumber, integrationNodesNumber - j));
+	marginalNodes[i][j] = 0.5 * (1.0 + p.x());
+	marginalWeights[i][j] = 0.5 * p.weight;
+      } // For j
   } // For i
   // Now, generate the nD rule over [0, 1]^n
   IndicesCollection allTuples(Tuples(discretization_).generate());