poly: Use a fixed number of DKW iterations for roots.

Also use a user supplied tolerance for detecting zeros.

Author: memmett

include/pfasst/quadrature/polynomial.hpp

@@ -33,7 +33,7 @@ namespace pfasst
       }
 
       Polynomial<CoeffT> normalize() const;
-      vector<CoeffT> roots() const;
+      vector<CoeffT> roots(size_t num_iterations=20, CoeffT ztol=1.0e-20) const;
       static Polynomial<CoeffT> legendre(const size_t order);
   };
 }  // ::pfasst

src/pfasst/quadrature/polynomial_impl.hpp

@@ -4,7 +4,6 @@
 #include <cassert>
 #include <cmath>
 #include <complex>
-#include <limits>
 using namespace std;
 
 
@@ -58,21 +57,21 @@ namespace pfasst
   }
 
   template<typename CoeffT>
-  vector<CoeffT> Polynomial<CoeffT>::roots() const
+  vector<CoeffT> Polynomial<CoeffT>::roots(size_t num_iterations, CoeffT ztol) const
   {
     assert(c.size() >= 1);
     size_t n = c.size() - 1;
 
     // initial guess
-    Polynomial<complex<CoeffT>> z0(n), z1(n);
+    vector<complex<CoeffT>> z0(n), z1(n);
     for (size_t j = 0; j < n; j++) {
       z0[j] = pow(complex<double>(0.4, 0.9), j);
       z1[j] = z0[j];
     }
 
     // durand-kerner-weierstrass iterations
-    Polynomial<CoeffT> p = normalize();
-    for (size_t k = 0; k < 100; k++) {
+    Polynomial<CoeffT> p = this->normalize();
+    for (size_t k = 0; k < num_iterations; k++) {
       complex<CoeffT> num, den;
       for (size_t i = 0; i < n; i++) {
         num = p.evaluate(z0[i]);
@@ -83,18 +82,12 @@ namespace pfasst
         }
         z0[i] = z0[i] - num / den;
       }
-
-      // converged?
-      CoeffT acc = 0.0;
-      for (size_t j = 0; j < n; j++) { acc += abs(z0[j] - z1[j]); }
-      if (acc < 2 * numeric_limits<CoeffT>::epsilon()) { break; }
-
       z1 = z0;
     }
 
     vector<CoeffT> roots(n);
     for (size_t j = 0; j < n; j++) {
-      roots[j] = (abs(z0[j]) < 4 * numeric_limits<CoeffT>::epsilon()) ? 0.0 : real(z0[j]);
+      roots[j] = abs(z0[j]) < ztol ? 0.0 : real(z0[j]);
     }
 
     sort(roots.begin(), roots.end());