FentonWave performance

[mgrouch]

diff --git a/FentonWave.h b/FentonWave.h
index abcdef1..1234567 100644
--- a/FentonWave.h
+++ b/FentonWave.h
@@ class FentonWave {
+    struct ResidualCache {
+        Eigen::Array<Real, N + 1, 1> eta;
+        Eigen::Array<Real, N + 1, 1> x_m;
+        std::array<Eigen::Array<Real, N, 1>, N + 1> kj_eta;
+        std::array<Eigen::Array<Real, N, 1>, N + 1> sinh_by_cosh_vals;
+        std::array<Eigen::Array<Real, N, 1>, N + 1> cosh_by_cosh_vals;
+        std::array<Eigen::Array<Real, N, 1>, N + 1> sin_jx;
+        std::array<Eigen::Array<Real, N, 1>, N + 1> cos_jx;
+    };

@@ void optimize(...) {
-        for (int iter = 0; iter < 500 && error > Real(1e-8); ++iter) {
-            Eigen::Matrix<Real, NU, 1> f = compute_residual(coeffs, H, k, D);
+        ResidualCache cache;
+        for (int iter = 0; iter < 500 && error > Real(1e-8); ++iter) {
+            Eigen::Matrix<Real, NU, 1> f = compute_residual(coeffs, H, k, D, cache);
             error = f.cwiseAbs().maxCoeff();
-            Eigen::Matrix<Real, NU, NU> J = compute_jacobian(coeffs, H, k, D);
+            Eigen::Matrix<Real, NU, NU> J = compute_jacobian(coeffs, H, k, D, cache);
@@ Eigen::Matrix<Real, StateDim, 1>
-    compute_residual(const Eigen::Matrix<Real, StateDim, 1>& coeffs, Real H, Real k, Real D) {
+    compute_residual(const Eigen::Matrix<Real, StateDim, 1>& coeffs, Real H, Real k, Real D, ResidualCache& cache) {
...
-        for (int m = 0; m <= N; ++m) {
-            Real x_m = Real(M_PI) * m / N;
-            Real eta_m = eta(m);
+        for (int m = 0; m <= N; ++m) {
+            Real x_m = Real(M_PI) * m / N;
+            Real eta_m = eta(m);
+            cache.eta(m) = eta_m;
+            cache.x_m(m) = x_m;
+            for (int j = 0; j < N; ++j) {
+                Real kj_eta = kj(j) * eta_m;
+                Real kj_D   = kj(j) * D;
+                cache.kj_eta[m](j) = kj_eta;
+                cache.sinh_by_cosh_vals[m](j) = sinh_by_cosh(kj_eta, kj_D);
+                cache.cosh_by_cosh_vals[m](j) = cosh_by_cosh(kj_eta, kj_D);
+                Real jx = j_cache(j) * x_m;
+                cache.sin_jx[m](j) = std::sin(jx);
+                cache.cos_jx[m](j) = std::cos(jx);
+            }
+        }

         for (int m = 0; m <= N; ++m) {
-            Real x_m = Real(M_PI) * m / N;
-            Real eta_m = eta(m);
-            Eigen::Array<Real, N, 1> kj_eta = kj * eta_m;
-            Eigen::Array<Real, N, 1> kj_D   = kj * D;
-            Eigen::Array<Real, N, 1> sin_jx = (j * x_m).sin();
-            Eigen::Array<Real, N, 1> cos_jx = (j * x_m).cos();
-            Eigen::Array<Real, N, 1> S1 = kj_eta.binaryExpr(kj_D, sinh_by_cosh);
-            Eigen::Array<Real, N, 1> C1 = kj_eta.binaryExpr(kj_D, cosh_by_cosh);
+            Real x_m = cache.x_m(m);
+            Real eta_m = cache.eta(m);
+            const auto& S1 = cache.sinh_by_cosh_vals[m];
+            const auto& C1 = cache.cosh_by_cosh_vals[m];
+            const auto& sin_jx = cache.sin_jx[m];
+            const auto& cos_jx = cache.cos_jx[m];
@@ Eigen::Matrix<Real, StateDim, StateDim>
-    compute_jacobian(const Eigen::Matrix<Real, StateDim, 1>& coeffs, Real H, Real k, Real D) {
+    compute_jacobian(const Eigen::Matrix<Real, StateDim, 1>& coeffs, Real H, Real k, Real D, const ResidualCache& cache) {
...
         for (int m = 0; m <= N; ++m) {
-            Real eta_m = eta(m);
-            Real x_m = Real(M_PI) * m / N;
+            Real eta_m = cache.eta(m);
+            Real x_m   = cache.x_m(m);
-            for (int j = 0; j < N; ++j) {
-                Real kj_eta = kj(j) * eta_m;
-                Real kj_D   = kj(j) * D;
-                S1(j) = sinh_by_cosh(kj_eta, kj_D);
-                C1(j) = cosh_by_cosh(kj_eta, kj_D);
-                Real jx = j_arr(j) * x_m;
-                S2(j) = std::sin(jx);
-                C2(j) = std::cos(jx);
-            }
+            S1 = cache.sinh_by_cosh_vals[m];
+            C1 = cache.cosh_by_cosh_vals[m];
+            S2 = cache.sin_jx[m];
+            C2 = cache.cos_jx[m];

[mgrouch]

--- FentonWave.h
+++ FentonWave.h
@@ using VelocityTerms = Eigen::Array<Real, N, 1>;
+    // ------------------------------------------------------------------------
+    // Per-iteration cache for residual & Jacobian
+    struct ResidualCache {
+      Eigen::Array<Real, N+1, 1>            x_m;
+      Eigen::Array<Real, N,   1>            kj, kjD;
+      std::array<Eigen::Array<Real, N,1>, N+1> sin_jx, cos_jx,
+                                               sinh_by_cosh_vals, cosh_by_cosh_vals;
+
+      ResidualCache(Real k, Real D, const Eigen::Array<Real,N,1>& j_arr) {
+        for (int m=0; m<=N; ++m) x_m(m) = Real(M_PI)*m/N;
+        kj  = j_arr * k;  kjD = kj * D;
+        for (int m=0; m<=N; ++m) for(int j=0;j<N;++j){
+          Real a = j_arr(j)*x_m(m);
+          sin_jx[m](j)=std::sin(a);
+          cos_jx[m](j)=std::cos(a);
+          sinh_by_cosh_vals[m](j)=0;
+          cosh_by_cosh_vals[m](j)=0;
+        }
+      }
+      void update_hyperbolic(const Eigen::Array<Real,N+1,1>& eta){
+        for(int m=0;m<=N;++m) for(int j=0;j<N;++j){
+          Real A=kj(j)*eta(m), B=kjD(j);
+          sinh_by_cosh_vals[m](j)=sinh_by_cosh(A,B);
+          cosh_by_cosh_vals[m](j)=cosh_by_cosh(A,B);
+        }
+      }
+    };
 
@@
-    void optimize(VectorF& B, Real& Q, Real& R, VectorF& eta, Real H, Real k, Real D) {
+    void optimize(VectorF& B, Real& Q, Real& R, VectorF& eta, Real H, Real k, Real D) {
@@
-      Real error = std::numeric_limits<Real>::max();
-      for (int iter = 0; iter < 500 && error > Real(1e-8); ++iter) {
+      Real error = std::numeric_limits<Real>::max();
+      ResidualCache cache{k,D,j_cache.array()};
+      for (int iter = 0; iter < 500 && error > Real(1e-8); ++iter) {
+        cache.update_hyperbolic(eta);
@@
-        Eigen::Matrix<Real, NU, 1> f = compute_residual(coeffs, H, k, D);
+        Eigen::Matrix<Real, NU, 1> f = compute_residual(coeffs, H, k, D, cache);
@@
-        Eigen::Matrix<Real, NU, NU> J = compute_jacobian(coeffs, H, k, D);
+        Eigen::Matrix<Real, NU, NU> J = compute_jacobian(coeffs, H, k, D, cache);
 
@@
-    Eigen::Matrix<Real, StateDim, 1>
-    compute_residual(const Eigen::Matrix<Real, StateDim, 1>& coeffs, Real H, Real k, Real D) {
+    Eigen::Matrix<Real, StateDim, 1>
+    compute_residual(const Eigen::Matrix<Real, StateDim, 1>& coeffs,
+                     Real H, Real k, Real D,
+                     const ResidualCache& cache) {
@@
-      for (int m = 0; m <= N; ++m) {
+      for (int m = 0; m <= N; ++m) {
@@
-        Eigen::Array<Real, N, 1> kj_eta = kj * eta_m;
-        Eigen::Array<Real, N, 1> kj_D   = kj * D;
-        Eigen::Array<Real, N, 1> S1     = kj_eta.binaryExpr(kj_D, sinh_by_cosh);
-        Eigen::Array<Real, N, 1> C1     = kj_eta.binaryExpr(kj_D, cosh_by_cosh);
+        auto const& S1 = cache.sinh_by_cosh_vals[m];
+        auto const& C1 = cache.cosh_by_cosh_vals[m];
+        auto const& sj = cache.sin_jx[m];
+        auto const& cj = cache.cos_jx[m];
 
@@
-    Eigen::Matrix<Real, StateDim, StateDim>
-    compute_jacobian(const Eigen::Matrix<Real, StateDim, 1>& coeffs, Real H, Real k, Real D) {
+    Eigen::Matrix<Real, StateDim, StateDim>
+    compute_jacobian(const Eigen::Matrix<Real, StateDim, 1>& coeffs,
+                     Real H, Real k, Real D,
+                     const ResidualCache& cache) {
@@
-      for (int m = 0; m <= N; ++m) {
+      for (int m = 0; m <= N; ++m) {
@@
-        for (int j = 0; j < N; ++j) {
-          S1(j) = sinh_by_cosh(kj(j)*eta_m, kj(j)*D);
-          C1(j) = cosh_by_cosh(kj(j)*eta_m, kj(j)*D);
-          S2(j) = std::sin(j_arr(j)*x_m);
-          C2(j) = std::cos(j_arr(j)*x_m);
-        }
+        auto const& S1 = cache.sinh_by_cosh_vals[m];
+        auto const& C1 = cache.cosh_by_cosh_vals[m];
+        auto const& S2 = cache.sin_jx[m];
+        auto const& C2 = cache.cos_jx[m];

[mgrouch]

Below are the exact hunks you need to apply to FentonWave.h. These show only the lines to remove (prefixed -) and the lines to add (prefixed +), with explicit types (no auto). Apply them in order.

⸻

1. In optimize(), build/update cache and change calls

--- FentonWave.h
+++ FentonWave.h
@@ void optimize(VectorF& B, Real& Q, Real& R, VectorF& eta, Real H, Real k, Real D) {
-      Real error = std::numeric_limits<Real>::max();
+      Real error = std::numeric_limits<Real>::max();
+      // build cache once per optimize() call
+      ResidualCache cache{k, D, j_cache.array()};
@@ for (int iter = 0; iter < 500 && error > Real(1e-8); ++iter) {
-        Eigen::Matrix<Real, NU, 1> f = compute_residual(coeffs, H, k, D);
+        // update cached sinh/cosh for current eta
+        cache.update(eta);
+        Eigen::Matrix<Real, NU, 1> f = compute_residual(coeffs, H, k, D, cache);
@@
-        Eigen::Matrix<Real, NU, NU> J = compute_jacobian(coeffs, H, k, D);
+        Eigen::Matrix<Real, NU, NU> J = compute_jacobian(coeffs, H, k, D, cache);

⸻

2. Change compute_residual signature

--- FentonWave.h
+++ FentonWave.h
@@
-    BigVector
-    compute_residual(const Eigen::Matrix<Real, StateDim, 1>& coeffs,
-                     Real H, Real k, Real D)
+    BigVector
+    compute_residual(const Eigen::Matrix<Real, StateDim, 1>& coeffs,
+                     Real H, Real k, Real D,
+                     const ResidualCache& cache)

⸻

3. Inside compute_residual, use cached arrays (replace inner loop setup)

Find the loop beginning for (int m = 0; m <= N; ++m) and remove the lines that explicitly recompute x_m, sin_jx, cos_jx, S1, C1. Insert the explicit cache usage instead:

--- FentonWave.h
+++ FentonWave.h
@@ compute_residual(...) {
-        Real x_m = Real(M_PI) * m / N;
-        Real eta_m = eta(m);
-
-        Eigen::Array<Real, N, 1> kj_eta = kj * eta_m;
-        Eigen::Array<Real, N, 1> kj_D   = kj * D;
-
-        Eigen::Array<Real, N, 1> sin_jx = (j * x_m).sin();
-        Eigen::Array<Real, N, 1> cos_jx = (j * x_m).cos();
-
-        Eigen::Array<Real, N, 1> S1 = kj_eta.binaryExpr(kj_D, [](Real a, Real b) {
-          return sinh_by_cosh(a, b);
-        });
-        Eigen::Array<Real, N, 1> C1 = kj_eta.binaryExpr(kj_D, [](Real a, Real b) {
-          return cosh_by_cosh(a, b);
-        });
+        // use cached trig & hyperbolic for this m
+        const Eigen::Array<Real, N, 1>& sj = cache.sin_jx[m];
+        const Eigen::Array<Real, N, 1>& cj = cache.cos_jx[m];
+        const Eigen::Array<Real, N, 1>& S1 = cache.sinh_vals[m];
+        const Eigen::Array<Real, N, 1>& C1 = cache.cosh_vals[m];
+
+        VelocityTerms SC = S1 * cj;
+        VelocityTerms CC = C1 * cj;
+        VelocityTerms SS = S1 * sj;

⸻

4) Change compute_jacobian signature

--- FentonWave.h
+++ FentonWave.h
@@
-    BigMatrix
-    compute_jacobian(const Eigen::Matrix<Real, StateDim, 1>& coeffs,
-                     Real H, Real k, Real D)
+    BigMatrix
+    compute_jacobian(const Eigen::Matrix<Real, StateDim, 1>& coeffs,
+                     Real H, Real k, Real D,
+                     const ResidualCache& cache)

⸻

5. Inside compute_jacobian, use cached arrays (replace inner loop setup)

Find the for (int m = 0; m <= N; ++m) in compute_jacobian and remove the explicit trig/hyperbolic loop. Insert:

--- FentonWave.h
+++ FentonWave.h
@@ compute_jacobian(...) {
-        Real eta_m = eta(m);
-        Real x_m   = Real(M_PI) * m / N;
-        Eigen::Array<Real, N, 1> S1, C1, S2, C2;
-
-        for (int j = 0; j < N; ++j) {
-          Real kj_eta = kj(j) * eta_m;
-          Real kj_D   = kj(j) * D;
-          S1(j) = sinh_by_cosh(kj_eta, kj_D);
-          C1(j) = cosh_by_cosh(kj_eta, kj_D);
-          Real jx    = j_arr(j) * x_m;
-          S2(j) = std::sin(jx);
-          C2(j) = std::cos(jx);
-        }
-
-        SC = S1 * C2;
-        SS = S1 * S2;
-        CC = C1 * C2;
-        CS = C1 * S2;
+        // use cached trig & hyperbolic for this m
+        const Eigen::Array<Real, N, 1>& sj = cache.sin_jx[m];
+        const Eigen::Array<Real, N, 1>& cj = cache.cos_jx[m];
+        const Eigen::Array<Real, N, 1>& S1 = cache.sinh_vals[m];
+        const Eigen::Array<Real, N, 1>& C1 = cache.cosh_vals[m];
+
+        VelocityTerms SC = S1 * cj;
+        VelocityTerms CC = C1 * cj;
+        VelocityTerms SS = S1 * sj;
+        VelocityTerms CS = C1 * sj;

⸻

Once you’ve made all five of these edits, rebuild. You will have:
• A single ResidualCache shared between residual & Jacobian.
• No auto—all types are explicit.
• Use of your VelocityTerms alias for SC, CC, SS, CS.
• No redundant inner‐loop trig or hyperbolic calls.

[mgrouch]

Here are the main places where you’re recomputing things that depend only on the fixed wave parameters (N, length, depth, height) and so can be lifted into the constructor (or a small “precompute” helper) and cached:
1. Harmonic index array

RealArray create_harmonic_indices() const {
return RealArray::LinSpaced(N + 1, 0, N);
}

You call create_harmonic_indices() in almost every surface‐derivative method.
Cache once as a member RealArray j0_N;.

2.	Denominators for the velocity‐term hyperbolics

const auto denom = (kj * depth).cosh();

In compute_velocity_terms you do

const auto kj = kj_cache.array();
const auto denom = (kj * depth).cosh();

but kj_cache and depth never change after construction.
Cache both

VectorN kj_cache; // already cached
VelocityTerms cosh_kjD_cache; // = (kj_cachedepth).cosh()
VelocityTerms sinh_kjD_cache; // = (kj_cachedepth).sinh()

3.	Grid of collocation points

VectorF x_nd = VectorF::LinSpaced(N + 1, 0, lam/2);

You use x_nd only in compute(). If you ever re‐solve for different wave heights but same geometry, cache x_nd and its cosines:

VectorF x_nd;
MatrixNxP cos_j_xnd; // cos(j * x_nd) for j=1…N

4.	Static DCT matrices and weights

FentonFFT::cosine_matrix();
FentonFFT::weights();

These are already static const, so you’re good here—don’t recompute them.

5.	Arrays of j⋅π/N for residual/jacobian loops

Inside both compute_residual and compute_jacobian you do

Real x_m = M_PI * m / N;
sin_jx = (j * x_m).sin();
cos_jx = (j * x_m).cos();

but the table of sin(j⋅π⋅m/N) and cos(j⋅π⋅m/N) for m=0…N is independent of the Newton steps.
Cache two (N+1)×N tables:

MatrixNxP sin_jm; // sin(j * mπ/N)
MatrixNxP cos_jm; // cos(j * mπ/N)

6.	Linear j array in residual

Array<Real,N,1> j = Array<Real,N,1>::LinSpaced(N,1, N);

You already have j_cache; use that everywhere instead of rebuilding j.

⸻

Putting it all together

In your constructor (or a small precompute() called once), you could add:

// 1. flat harmonic indices 0…N
j0_N = RealArray::LinSpaced(N+1, 0, N);

// 2. depth‐dependent hyperbolic denominators
{
auto kj = kj_cache.array();
cosh_kjD_cache = (kj * depth).cosh();
sinh_kjD_cache = (kj * depth).sinh();
}

// 3. collocation grid and its cosines
x_nd = VectorF::LinSpaced(N+1, 0, lam/2);
for(int j=1; j<=N; ++j)
cos_j_xnd.row(j-1) = (x_nd * Real(j)).array().cos();

// 5. sin/cos tables for residual/jacobian
for(int m=0; m<=N; ++m){
Real xm = Real(M_PI)*m/N;
sin_jm.row(m) = (j_cache.array() * xm).sin();
cos_jm.row(m) = (j_cache.array() * xm).cos();
}

Then in your hot loops:
• Replace create_harmonic_indices() with j0_N.
• Replace (kjdepth).cosh() / .sinh() with the cached arrays.
• Replace on‐the‐fly sin(jx_m) / cos(j*x_m) with table lookups in sin_jm[m]/cos_jm[m].
• Use j_cache instead of rebuilding j = LinSpaced(1…N).

This should eliminate all the per‐call allocations, transforms, and heavy transcendental calls that only depend on the geometry and not on the Newton‐iteration state.