Optimize GridInterp tensor-product cubic interpolation

[luwang00]

Feature or improvement description

Optimizes the tensor-product cubic interpolation routines in the NWTC Library's GridInterp module for improved runtime performance, while preserving the public interface and functionality. The results are numerically equivalent to the original — differences arise only at floating-point truncation level due to summation regrouping, not from any change in the interpolation math.

Key optimizations:

• Setup routines (3D/4D/N): weight tensors are built with array syntax over the contiguous first dimension so the inner writes vectorize, and shared subexpressions are hoisted (e.g. Nkl = N1D(k,3)*N1D(l,4)).
• Scalar & vector interpolation: the reduction over the innermost dimension accumulates into four independent partial sums, breaking the serial add dependency chain (latency hiding) and exposing instruction-level parallelism; index lookups are hoisted out of the inner loop.
• Vector routines: the component (vi) loop is moved outermost for unit-stride access with per-component accumulation.
• Slope routines (N/S): both derivatives are accumulated with dual accumulators, reusing a single data load per grid point instead of loading the same value twice.

Rationale:
gfortran does not reassociate reduction sums by default, so the original serial add chain was latency-bound. The multi-accumulator regrouping is the primary win.

Related issue, if one exists

N/A

Impacted areas of the software

GridInterp.f90 (only file changed)
Downstream consumers exercised: SeaState (SeaSt_WaveField) and HydroDyn/WAMIT (WAMIT_Interp)

Additional supporting information

The amortization design of the setup routines (precomputing full weight tensors reused across many field interpolations at the same position) was preserved intentionally. No interface, argument, or return-type changes; the diff is 202 insertions / 112 deletions in a single file.

Generative AI usage

The optimization was developed with AI assistance (analysis, refactoring, and validation).

Co-authored-by: Microsoft Copilot [copilot@microsoft.com]
Co-authored-by: Anthropic Claude [claude@anthropic.com]

Test results, if applicable
No change to existing r-test expected, pending regression test completion.

[Copilot]

Pull request overview

This PR optimizes the NWTC Library GridInterp tensor-product cubic interpolation implementation to improve runtime performance while keeping the public API unchanged, primarily by increasing vectorization opportunities and reducing reduction-chain latency in inner loops.

Changes:

• Rewrites 3D/4D/N setup routines to build weight tensors using whole-slice array syntax and hoisted shared subexpressions.
• Refactors scalar and vector interpolation kernels to use multiple accumulators and hoist index lookups out of inner loops.
• Refactors slope (S/N) routines to reuse a single grid-point load for both derivative accumulations.

Suppressed comments (3)

modules/nwtc-library/src/GridInterp.f90:411

• In GridInterpSetupN, the local variable i is declared but no longer used after replacing the i loop with whole-slice assignments into m%N4D. Dropping it will avoid unused-variable warnings.

   ! Need two sets of weights for d(.)/dx and d(.)/dy. Borrow m%N4D for this.
   do k = 1,4
      do j = 1,4
         m%N4D(:,j,k,1) = N1D(:,1) * (N1Ddx(j,2)*N1D  (k,3))
         m%N4D(:,j,k,2) = N1D(:,1) * (N1D  (j,2)*N1Ddx(k,3))

modules/nwtc-library/src/GridInterp.f90:841

• GridInterp4DVecNR8 no longer initializes the result array before the vi loop. If vDim is 0, the loop won't run and the function result is undefined, whereas the previous implementation defined it to 0.0 via whole-array initialization.

   ! interpolate
   do vi = 1,vDim
      acc = 0.0_DbKi

modules/nwtc-library/src/GridInterp.f90:365

• i is declared in GridInterpSetup4D but no longer used after converting the inner i loop to an array slice assignment. Removing it avoids unused-variable warnings.

      do k = 1,4
         Nkl = N1D(k,3)*N1D(l,4)
         do j = 1,4
            m%N4D(:,j,k,l) = N1D(:,1) * (N1D(j,2)*Nkl)

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.