Use a callable template param in DmpCeff

CMakeLists.txt

@@ -73,7 +73,6 @@ set(STA_SOURCE
   dcalc/DelayCalcBase.cc
   dcalc/DmpCeff.cc
   dcalc/DmpDelayCalc.cc
-  dcalc/FindRoot.cc
   dcalc/GraphDelayCalc.cc
   dcalc/LumpedCapDelayCalc.cc
   dcalc/NetCaps.cc

dcalc/DmpCeff.cc

@@ -92,13 +92,14 @@ gateModelRd(const LibertyCell *cell,
 	    double c1,
 	    const Pvt *pvt,
 	    bool pocv_enabled);
+template <typename Eval>
 static void
 newtonRaphson(const int max_iter,
 	      double x[],
 	      const int n,
 	      const double x_tol,
 	      // eval(state) is called to fill fvec and fjac.
-	      function<void ()> eval,
+	      const Eval &eval,
 	      // Temporaries supplied by caller.
 	      double *fvec,
 	      double **fjac,
@@ -494,9 +495,9 @@ DmpAlg::findVoCrossing(double vth,
                        double t_lower,
                        double t_upper)
 {
-  FindRootFunc vo_func = [&] (double t,
-                              double &y,
-                              double &dy) {
+  const auto vo_func = [&] (double t,
+                            double &y,
+                            double &dy) {
     double vo, vo_dt;
     Vo(t, vo, vo_dt);
     y = vo - vth;
@@ -612,9 +613,9 @@ DmpAlg::findVlCrossing(double vth,
                        double t_lower,
                        double t_upper)
 {
-  FindRootFunc vl_func = [&] (double t,
-                              double &y,
-                              double &dy) {
+  const auto vl_func = [&] (double t,
+                            double &y,
+                            double &dy) {
     double vl, vl_dt;
     Vl(t, vl, vl_dt);
     y = vl - vth;
@@ -1278,19 +1279,22 @@ DmpZeroC2::voCrossingUpperBound()
 // x_tol is percentage that all changes in x must be less than (1.0 = 100%).
 // Eval(state) is called to fill fvec and fjac (returns false if fails).
 // Return error msg on failure.
+// Eval should be callable like void(void)
+template <typename Eval>
 static void
 newtonRaphson(const int max_iter,
 	      double x[],
 	      const int size,
 	      const double x_tol,
-	      function<void ()> eval,
+	      const Eval &eval,
 	      // Temporaries supplied by caller.
 	      double *fvec,
 	      double **fjac,
 	      int *index,
 	      double *p,
 	      double *scale)
 {
+  static_assert(std::is_invocable_v<Eval>, "Must be a callable with no arguments");
   for (int k = 0; k < max_iter; k++) {
     eval();
     for (int i = 0; i < size; i++)

dcalc/FindRoot.hh

@@ -16,33 +16,95 @@
 
 #pragma once
 
-#include <functional>
+#include <algorithm> // abs
 
 namespace sta {
 
-typedef const std::function<void (double x,
-                                  // Return values.
-                                  double &y,
-                                  double &dy)> FindRootFunc;
+using std::abs;
 
+// FindRootFunc should be callable like void(double, double&, double&)
+template <typename FindRootFunc>
 double
-findRoot(FindRootFunc func,
+findRoot(const FindRootFunc &func,
 	 double x1,
 	 double x2,
 	 double x_tol,
 	 int max_iter,
          // Return value.
-         bool &fail);
+         bool &fail)
+{
+  double y1, y2, dy1;
+  func(x1, y1, dy1);
+  func(x2, y2, dy1);
+  return findRoot(func, x1, y1, x2, y2, x_tol, max_iter, fail);
+}
 
+// FindRootFunc should be callable like void(double, double&, double&)
+template <typename FindRootFunc>
 double
-findRoot(FindRootFunc func,
+findRoot(const FindRootFunc &func,
 	 double x1,
 	 double y1,
          double x2,
 	 double y2,
 	 double x_tol,
 	 int max_iter,
          // Return value.
-         bool &fail);
+         bool &fail)
+{
+  if ((y1 > 0.0 && y2 > 0.0) || (y1 < 0.0 && y2 < 0.0)) {
+    // Initial bounds do not surround a root.
+    fail = true;
+    return 0.0;
+  }
+
+  if (y1 == 0.0) {
+    fail = false;
+    return x1;
+  }
+
+  if (y2 == 0.0) {
+    fail = false;
+    return x2;
+  }
+
+  if (y1 > 0.0)
+    // Swap x1/x2 so func(x1) < 0.
+    std::swap(x1, x2);
+  double root = (x1 + x2) * 0.5;
+  double dx_prev = abs(x2 - x1);
+  double dx = dx_prev;
+  double y, dy;
+  func(root, y, dy);
+  for (int iter = 0; iter < max_iter; iter++) {
+    // Newton/raphson out of range.
+    if ((((root - x2) * dy - y) * ((root - x1) * dy - y) > 0.0)
+	// Not decreasing fast enough.
+	|| (abs(2.0 * y) > abs(dx_prev * dy))) {
+      // Bisect x1/x2 interval.
+      dx_prev = dx;
+      dx = (x2 - x1) * 0.5;
+      root = x1 + dx;
+    }
+    else {
+      dx_prev = dx;
+      dx = y / dy;
+      root -= dx;
+    }
+    if (abs(dx) <= x_tol * abs(root)) {
+      // Converged.
+      fail = false;
+      return root;
+    }
+
+    func(root, y, dy);
+    if (y < 0.0)
+      x1 = root;
+    else
+      x2 = root;
+  }
+  fail = true;
+  return root;
+}
 
 } // namespace