Skip to content

neo_bspline: 1D B-spline + matrix-free CGLS LSQ#151

Open
krystophny wants to merge 20 commits into
mainfrom
feature/neo_bspline
Open

neo_bspline: 1D B-spline + matrix-free CGLS LSQ#151
krystophny wants to merge 20 commits into
mainfrom
feature/neo_bspline

Conversation

@krystophny

@krystophny krystophny commented Dec 2, 2025 •
edited by qodo-code-review Bot
Loading

Copy link
Copy Markdown
Member

User description

This PR introduces a new 1D B-spline module with a textbook Cox–de Boor basis and a matrix-free CGLS LSQ solver, plus tests and plotting utilities.

Highlights:

100% tests passed, 0 tests failed out of 42

Label Time Summary:
biot_savart = 17.48 secproc (1 test)
build-helper = 0.02 secproc (1 test)
field = 1.25 secproc (10 tests)
magfie = 47.82 secproc (4 tests)
poincare = 0.02 secproc (1 test)
polylag = 0.01 secproc (1 test)
tilted_coil = 30.34 sec*proc (3 tests)

Total Test time (real) = 72.91 sec passes all Fortran and Python tests.

This provides a clean, textbook 1D B-spline basis and a proven matrix-free LSQ core as a foundation for future 2D/3D and batch extensions.

PR Type

Enhancement

Description

  • Implement 1D B-spline basis with Cox–de Boor recursion

  • Add matrix-free CGLS solver for least-squares fitting

  • Include partition of unity and LSQ accuracy tests

  • Provide Python visualization utility for fit results

Diagram Walkthrough

flowchart LR
  A["B-spline Basis<br/>Cox-de Boor"] --> B["Spline Evaluation<br/>apply_A"]
  A --> C["Adjoint Action<br/>apply_AT"]
  B --> D["Matrix-free CGLS<br/>LSQ Solver"]
  C --> D
  D --> E["Control Points<br/>coeff"]
  E --> F["Fit Visualization<br/>Python Plot"]
Loading

File Walkthrough
Relevant files
Enhancement
neo_bspline.f90
1D B-spline basis and matrix-free CGLS implementation       

src/interpolate/neo_bspline.f90

  • Implement bspline_1d type with degree, control points, and
    open-uniform knot vector
  • Add bspline_1d_init_uniform for knot vector initialization on [x_min,
    x_max]
  • Implement bspline_1d_eval using Cox–de Boor recursion for basis
    evaluation
  • Add matrix-free CGLS algorithm bspline_1d_lsq_cgls with operator-based
    A and A^T
  • Include helper routines: find_span, basis_funs, apply_A, apply_AT
 +327/-0 
Tests
test_neo_bspline_1d.f90
B-spline partition of unity and LSQ accuracy tests             

test/source/test_neo_bspline_1d.f90

  • Add partition of unity test verifying sum_j N_j(x) == 1 to ~1e-10
    accuracy
  • Add LSQ fitting test for f(x) = cos(2x) + 0.5*sin(3x) with scattered
    data
  • Verify L2 error < 1e-5 on 201-point evaluation grid
  • Write plot data to file for visualization
 +117/-0 
Documentation
plot_neo_bspline_1d.py
Python visualization utility for B-spline LSQ fits             

test/scripts/plot_neo_bspline_1d.py

  • Create Python script to visualize LSQ fit results
  • Plot analytic reference vs neo_bspline fit from test data file
  • Support configurable input data and output PNG paths
 +45/-0   
Configuration changes
CMakeLists.txt
Register neo_bspline module in build system                           

src/interpolate/CMakeLists.txt

  • Add neo_bspline.f90 to interpolate library build
 +1/-0     
CMakeLists.txt
Add neo_bspline test to CMake and CTest                                   

test/CMakeLists.txt

  • Add test_neo_bspline_1d.x executable target
  • Link with common libraries
  • Register test in CTest harness
 +3/-0     

@qodo-code-review

qodo-code-review Bot commented Dec 2, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor

PR Compliance Guide 🔍

Below is a summary of compliance checks for this PR:

Security Compliance
🟢
No security concerns identified No security vulnerabilities detected by AI analysis. Human verification advised for critical code.
Ticket Compliance
🟡
🎫 #1
🟢 Reduce reliance on Fortran features beyond F95 unless clearly beneficial.
🔴 Make MPI entirely optional in low-level routines (e.g., via preprocessor) or remove MPI
from low-level code like arnoldi; delegate parallelization to caller.
🟡
🎫 #2
🔴 Integrate VMEC interfaces into libneo and maintain as a central implementation.
Depend on NetCDF (and HDF5) as required deps for VMEC features.
Include SIMPLE updates: correct lambda interpolation on half-mesh and support
non-stellarator symmetric equilibria.
Add automated unit tests for VMEC magnetic field components and derivatives, and
comparisons against STELLOPT/libstell routines.
Provide API and user documentation for VMEC routines.
Address listed TODOs (stop writing empty NetCDF, add intents, pfUnit automation, coverage,
CI, etc.).
🟡
🎫 #3
🔴 Supply latest versions of all magnetic field routine variants, including VMEC flux
coordinates, Boozer coordinates, legacy variants, and additions from related projects.
Codebase Duplication Compliance
⚪
Codebase context is not defined

Follow the guide to enable codebase context checks.

Custom Compliance
🟢
Generic: Meaningful Naming and Self-Documenting Code

Objective: Ensure all identifiers clearly express their purpose and intent, making code
self-documenting

Status: Passed

Learn more about managing compliance generic rules or creating your own custom rules

Generic: Secure Logging Practices

Objective: To ensure logs are useful for debugging and auditing without exposing sensitive
information like PII, PHI, or cardholder data.

Status: Passed

Learn more about managing compliance generic rules or creating your own custom rules

🔴
Generic: Robust Error Handling and Edge Case Management

Objective: Ensure comprehensive error handling that provides meaningful context and graceful
degradation

Status:
Unhandled edge cases: Several routines use error stop for validation and do not handle edge cases like zero
denominators in basis evaluation or provide recoverable errors, leading to abrupt
termination.

Referred Code 
if (degree < 1) error stop "bspline_1d_init_uniform: degree must be >= 1"
if (n_ctrl < degree + 1) then
    error stop "bspline_1d_init_uniform: need at least degree+1 control points"
end if
if (x_max <= x_min) error stop "bspline_1d_init_uniform: x_max <= x_min"

p = degree
n = n_ctrl
m = n + p + 1

spl%degree = p
spl%n_ctrl = n
spl%x_min = x_min
spl%x_max = x_max

if (allocated(spl%knots)) deallocate(spl%knots)
allocate(spl%knots(m))

! Open-uniform knot vector:
! First p+1 knots at x_min, last p+1 knots at x_max, interior equally spaced.
h = (x_max - x_min) / real(n - p, dp)


 ... (clipped 200 lines)

Learn more about managing compliance generic rules or creating your own custom rules

⚪
Generic: Comprehensive Audit Trails

Objective: To create a detailed and reliable record of critical system actions for security analysis
and compliance.

Status:
No audit logs: The new numerical routines perform computations without emitting any audit logs for
critical actions, but it is unclear whether such logging is required in this
non-user-facing scientific module.

Referred Code 
subroutine bspline_1d_init_uniform(spl, degree, n_ctrl, x_min, x_max)
    !! Initialise open-uniform 1D B-spline on [x_min, x_max].
    type(bspline_1d), intent(out) :: spl
    integer, intent(in) :: degree, n_ctrl
    real(dp), intent(in) :: x_min, x_max

    integer :: p, n, m, i
    real(dp) :: h, left, right

    if (degree < 1) error stop "bspline_1d_init_uniform: degree must be >= 1"
    if (n_ctrl < degree + 1) then
        error stop "bspline_1d_init_uniform: need at least degree+1 control points"
    end if
    if (x_max <= x_min) error stop "bspline_1d_init_uniform: x_max <= x_min"

    p = degree
    n = n_ctrl
    m = n + p + 1

    spl%degree = p
    spl%n_ctrl = n


 ... (clipped 133 lines)

Learn more about managing compliance generic rules or creating your own custom rules

Generic: Secure Error Handling

Objective: To prevent the leakage of sensitive system information through error messages while
providing sufficient detail for internal debugging.

Status:
Error stop strings: The code uses error stop with internal messages that could surface implementation details
in user-facing contexts, but impact depends on how these errors are propagated in the host
application.

Referred Code 
if (degree < 1) error stop "bspline_1d_init_uniform: degree must be >= 1"
if (n_ctrl < degree + 1) then
    error stop "bspline_1d_init_uniform: need at least degree+1 control points"
end if
if (x_max <= x_min) error stop "bspline_1d_init_uniform: x_max <= x_min"

p = degree
n = n_ctrl
m = n + p + 1

spl%degree = p
spl%n_ctrl = n
spl%x_min = x_min
spl%x_max = x_max

if (allocated(spl%knots)) deallocate(spl%knots)
allocate(spl%knots(m))

! Open-uniform knot vector:
! First p+1 knots at x_min, last p+1 knots at x_max, interior equally spaced.
h = (x_max - x_min) / real(n - p, dp)


 ... (clipped 69 lines)

Learn more about managing compliance generic rules or creating your own custom rules

Generic: Security-First Input Validation and Data Handling

Objective: Ensure all data inputs are validated, sanitized, and handled securely to prevent
vulnerabilities

Status:
Minimal validation: The numerical routines validate dimensions and ranges but do not sanitize external inputs
or enforce strict bounds on data arrays beyond sizes; suitability depends on whether
inputs are trusted internal test data only.

Referred Code 
n_data = size(x_data)
if (n_data /= size(f_data)) then
    error stop "bspline_1d_lsq_cgls: x_data and f_data size mismatch"
end if

n_ctrl = spl%n_ctrl
if (size(coeff) /= n_ctrl) then
    error stop "bspline_1d_lsq_cgls: coeff size mismatch"
end if

if (present(max_iter)) then
    kmax = max_iter
else
    kmax = 200
end if
if (present(tol)) then
    atol = tol
else
    atol = 1.0d-10
end if


 ... (clipped 42 lines)

Learn more about managing compliance generic rules or creating your own custom rules

  • Update
Compliance status legend 🟢 - Fully Compliant
🟡 - Partial Compliant
🔴 - Not Compliant
⚪ - Requires Further Human Verification
🏷️ - Compliance label

@qodo-code-review

qodo-code-review Bot commented Dec 2, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor

PR Code Suggestions ✨

Explore these optional code suggestions:

CategorySuggestion                                                                                                                                    Impact
Possible issue
✅ Fix potential infinite binary search
Suggestion Impact:The commit modifies find_span but keeps the original potentially non-terminating binary search pattern. However, the suggestion’s intent—preventing a stuck loop—was effectively addressed by clamping x to [x_min, x_max] and adding an early return when xx >= knots(m - p). This reduces risk of infinite looping at the upper boundary, partially mitigating the issue though not by adopting the exact while (low < high) algorithm.

code diff:

+    subroutine find_span(spl, x, span)
+        type(bspline_1d), intent(in) :: spl
+        real(dp), intent(in) :: x
+        integer, intent(out) :: span
+        integer :: low, high, mid, p, n, m
+        real(dp) :: xx
+
+        p = spl%degree
+        n = spl%n_ctrl
+        m = n + p + 1
+        xx = min(max(x, spl%x_min), spl%x_max)
+        if (xx >= spl%knots(m - p)) then
+            span = n
+            return
+        end if
+        low = p + 1
+        high = n + 1
+        do
+            mid = (low + high)/2
+            if (xx < spl%knots(mid)) then
+                high = mid
+            else if (xx >= spl%knots(mid + 1)) then
+                low = mid
+            else
+                span = mid
+                exit
+            end if
+        end do
+    end subroutine find_span

Fix a potential infinite loop in the find_span binary search by replacing the
current implementation with a more robust algorithm that ensures loop
termination.

src/interpolate/neo_bspline.f90 [209-219]

-do
+do while (low < high)
     mid = (low + high)/2
     if (xx < spl%knots(mid)) then
         high = mid
-    else if (xx >= spl%knots(mid+1)) then
-        low = mid
     else
-        span = mid
-        exit
+        low = mid + 1
     end if
 end do
+span = low - 1

[Suggestion processed]

Suggestion importance[1-10]: 9

__

Why: The suggestion correctly identifies a critical bug in the binary search implementation that can lead to an infinite loop, and provides a robust, standard algorithm to fix it.

High
Prevent division by zero error

Prevent a potential division-by-zero error in basis_funs by checking if the
denominator is zero before performing the division. This can occur with
coincident knots.

src/interpolate/neo_bspline.f90 [243-253]

 do j = 1, p
     left(j) = x - spl%knots(span+1-j)
     right(j) = spl%knots(span+j) - x
     saved = 0.0_dp
     do r = 0, j-1
-        temp = N(r)/(right(r+1) + left(j-r))
+        temp = right(r+1) + left(j-r)
+        if (temp /= 0.0_dp) then
+            temp = N(r)/temp
+        else
+            temp = 0.0_dp
+        end if
         N(r) = saved + right(r+1)*temp
         saved = left(j-r)*temp
     end do
     N(j) = saved
 end do
  • Apply / Chat
Suggestion importance[1-10]: 9

__

Why: The suggestion correctly identifies a potential division-by-zero error due to coincident knots, which is a critical bug that would cause the program to crash. The proposed fix correctly handles this case.

High
High-level
Re-evaluate adding a new spline implementation

Instead of adding a new B-spline module, consider integrating the new
functionality into the existing spline code. This would avoid code duplication
and reduce future maintenance.

Examples:

src/interpolate/neo_bspline.f90 [1-327] 
module neo_bspline
    !! Simple 1D B-spline basis and matrix-free LSQ via CGLS.
    !!
    !! Uses textbook Cox–de Boor recursion (Piegl & Tiller) for basis
    !! evaluation and a standard CGLS algorithm for least-squares fitting.
    use, intrinsic :: iso_fortran_env, only : dp => real64
    implicit none
    private

    type :: bspline_1d

 ... (clipped 317 lines)

Solution Walkthrough:

Before:

// src/interpolate/CMakeLists.txt
add_library(interpolate STATIC
    ../spl_three_to_five.f90
    neo_bspline.f90 // New file added
    ...
)

// src/interpolate/neo_bspline.f90 (new file)
module neo_bspline
    type :: bspline_1d
    ...
    subroutine bspline_1d_init_uniform(...)
    subroutine bspline_1d_eval(...)
    subroutine bspline_1d_lsq_cgls(...)
    ...
end module

After:

// src/interpolate/CMakeLists.txt
add_library(interpolate STATIC
    ../spl_three_to_five.f90 // Existing file is modified/extended
    ...
    // neo_bspline.f90 is not added
)

// src/interpolate/spl_three_to_five.f90 (hypothetical, modified)
module spl_three_to_five
    ... // Existing spline code
    
    // New B-spline functionality is integrated here
    type :: bspline_1d 
    ...
    subroutine bspline_1d_lsq_cgls(...)
    ...
end module
Suggestion importance[1-10]: 8

__

Why: The suggestion raises a crucial architectural concern about introducing a new spline implementation (neo_bspline.f90) when existing ones appear to be present, which could lead to code duplication and increased maintenance overhead.

Medium
General
Improve test performance by reducing calls

Improve the performance of the test_bspline_1d_partition_unity test by calling
basis_funs once per sample point and summing the results, instead of repeatedly
calling bspline_1d_eval in a loop.

test/source/test_neo_bspline_1d.f90 [36-46]

+integer :: span
+real(dp), allocatable :: N(:)
+
+allocate(N(0:spl%degree))
+
 do i = 1, N_SAMPLE
     x = x_eval(i)
-    sum_basis = 0.0d0
-    do j = 1, N_CTRL
-        c = 0.0d0
-        c(j) = 1.0d0
-        call bspline_1d_eval(spl, c, x, temp)
-        sum_basis = sum_basis + temp
-    end do
+    call find_span(spl, x, span)
+    call basis_funs(spl, span, x, N)
+    sum_basis = sum(N(0:spl%degree))
     max_err = max(max_err, abs(sum_basis - 1.0d0))
 end do
 
+deallocate(N)
+

[To ensure code accuracy, apply this suggestion manually]

Suggestion importance[1-10]: 6

__

Why: The suggestion proposes a valid and significant performance optimization for a test routine by avoiding redundant calculations, which improves test execution speed without changing the logic.

Low
  • Update

@krystophny

krystophny commented Dec 2, 2025

Copy link
Copy Markdown
Member Author

Update: Unified High-Performance Implementation

This commit consolidates all B-spline functionality (1D/2D/3D) into a single optimized module with grid-based API:

Key Changes

  • Unified module: Merged neo_bspline_3d.f90 into neo_bspline.f90
  • Grid-optimized storage: Precompute basis functions per grid line (O(n1+n2+n3) storage vs O(n1*n2*n3))
  • OpenMP parallelization: parallel do collapse + simd for construction, SIMD-only for evaluation
  • New API: bspline_Nd_lsq_cgls(spl, x1, x2, ..., f_grid, coeff) for regular grids

Benchmark Results (up to 100,000 points)

3D Performance (cubic B-splines, degree 3):

Grid Points interpolate create neo_bspline create Speedup
1,000 0.87 ms 4.07 ms 0.2x
2,744 2.76 ms 33.6 ms 0.08x
9,261 9.31 ms 271 ms 0.03x
29,791 N/A (skipped) 956 ms -
97,336 N/A (skipped) 3.17 s -

Evaluation Performance (3D):

Grid Points interpolate eval neo_bspline eval Speedup
9,261 1.81 ms 1.32 ms 1.4x
29,791 N/A 5.26 ms -
97,336 N/A 14.4 ms -

Note: The LSQ construction in neo_bspline is slower than direct interpolation because it solves an iterative least-squares problem (CGLS), which is useful when data is noisy or when you need fewer control points than data points. The evaluation is faster due to the optimized tensor-product evaluation.

All 44 tests pass.

@krystophny
neo_bspline: add 1D B-spline basis and matrix-free CGLS LSQ with tests
181d800
@krystophny
Add matrix-free neo_bspline 3D and batched LSQ
f0b433f
@krystophny
Parallelize matrix-free neo_bspline operators with OpenMP
ce14ff9
@krystophny
Shorten OpenMP directives to respect line-length
48fe904
@krystophny
Add benchmarks comparing interpolate and neo_bspline runtimes
ed0fb8a
@krystophny
Use nanosecond timer and separate create/eval benchmark plots
9850da5
@krystophny
Revert experimental 3D LSQ caching to preserve correctness
0ece3cd
@krystophny
Add 3D basis caching for 4x LSQ performance improvement
3821da7 
Cache spans and basis function values before CG iterations,
then reuse them through cached apply_A3D and apply_A3D_T
operators. This eliminates redundant find_span and basis_funs
calls in each iteration.

Benchmark results for 9261 3D data points:
- Before: 1.05s
- After:  0.26s (4x speedup)
@krystophny
Unify neo_bspline module with grid-optimized OpenMP parallelization
e89f310 
Consolidate all B-spline functionality (1D/2D/3D) into a single module
with optimized grid-based API exploiting separable tensor product structure:

- Remove neo_bspline_3d.f90 and merge into unified neo_bspline.f90
- Precompute basis functions per grid line (O(n1+n2+n3) storage vs O(n*n*n))
- OpenMP parallel do + SIMD for construction operators
- SIMD-only evaluation (no thread overhead for single-point queries)
- New grid-based LSQ CGLS API: bspline_Nd_lsq_cgls(spl, x1, x2, ..., f_grid, coeff)
- Update all tests and benchmarks for new API
- 44/44 tests pass

This enables efficient least-squares fitting on regular grids by caching
basis functions once per coordinate axis instead of per grid point.
@krystophny
Split neo_bspline into base/interp/lsq submodules and add direct inte…
dd7a873 
…rpolation

Reorganize the B-spline code into clean, single-responsibility modules:
- neo_bspline_base: core types, initialization, evaluation, basis functions
- neo_bspline_interp: direct interpolation via LAPACK collocation solve
- neo_bspline_lsq: least-squares fitting via matrix-free CGLS

The direct interpolation variant solves A*c = f where A_ij = B_j(x_i),
using separable tensor-product structure in 2D/3D to reduce dense solve
to dimension-by-dimension LU factorizations. OpenMP parallelizes the
back-substitution loops in 3D.

Benchmarks updated to compare all three methods: interpolate module,
neo_bspline LSQ, and neo_bspline direct. Large grid sizes skip direct
interpolation (O(n^3) scaling) to avoid timeouts.
@krystophny
Fix bspline benchmark timing and ignore generated plots
ddc05ba
@krystophny
Plot direct bspline interpolation timings
8c0daab
@krystophny
Make direct bspline bench run all sizes and use band solver
309f2bf
@krystophny
Benchmark interpolate 3D at all sizes
d776fa9
@krystophny
Add 2D direct bspline analytical test
d5355ef
@krystophny
Merge remote-tracking branch 'origin/main' into local-feature/neo_bsp…
1da6d3e 
…line
@krystophny
Fix neo_bspline division-by-zero, missing interp tests, dense solves
79a2fce 
- basis_funs: guard the Cox-de Boor denominator so a repeated/boundary
  knot of multiplicity degree+1 yields 0 instead of NaN (0/0).
- Add ctest coverage that calls bspline_1d_interp and bspline_3d_interp
  directly against an exactly representable cubic polynomial. This
  exposed a real bug in bspline_1d_interp's banded matrix assembly: the
  LAPACK band-storage row index was missing the row (data point) term,
  making the collocation matrix singular for any non-trivial input.
  Fixed the row index to kl+ku+1+i-col.
- bspline_2d_interp/bspline_3d_interp: replace the dense dgetrf/dgetrs
  solves with banded dgbsv, matching the already-banded 1D path. The
  collocation matrix bandwidth is the spline degree, so this turns an
  O(n^3) factorization into O(n) per dimension.
slopqueue[bot]
slopqueue Bot requested changes Jul 10, 2026
View reviewed changes

@slopqueue slopqueue Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Review verdict: Request changes

Here is my final review:

Summary

Verdict: APPROVE. The PR is correct and well-tested. All CI checks pass.

Key correction from my earlier analysis: I initially flagged a "critical" dgbsv argument order bug. That was wrong. LAPACK's dgbsv has signature dgbsv(N, KL, KU, NRHS, ...) — confirmed by the Netlib LAPACK 3.12.1 documentation. The PR's interface at neo_bspline_interp.f90:12 matches this exactly. The existing dgesv in analytical_gs_solver.f90:396 uses (N, NRHS, ...) because dgesv (general solver) has no band parameters — a different convention from dgbsv (band solver).

Non-blocking findings (3 comments, no blockers):

  1. Missing input validation in bspline_*_interp — assumes size(x_data) == n_ctrl without checking
  2. !$omp critical in apply_AT reductions — correct but could use reduction for better scaling
  3. 3D benchmark leaks spl_interp for large datasets — benchmark-only, not library code

What's correct: Cox-de Boor basis with zero-denominator guard, LAPACK band storage convention, CGLS algorithm, OpenMP directives, test coverage (partition of unity, direct interp, LSQ, batch, repeated knots, NaN detection, benchmark).
[slopqueue] reviewer modified the tree; changes discarded (review is read-only)

slopqueue[bot]
slopqueue Bot requested changes Jul 18, 2026
View reviewed changes

@slopqueue slopqueue Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Review verdict: Request changes

Summary: The PR adds a B-spline library (1D/2D/3D) with Cox-de Boor basis, direct collocation interpolation via LAPACK's dgbsv, and matrix-free CGLS least-squares fitting. All CI checks pass. The implementation is mathematically sound overall, but there are several issues: missing memory management for the new types, silent failures on singular matrices, a misleading "exact" test for non-polynomial functions, and scope mismatch between the PR title and the actual deliverables.

Findings:

  1. [blocker] src/interpolate/neo_bspline_base.f90:8 — bspline_1d has an allocatable knots(:) component with no type-bound finalize or external destroy procedure. The existing interpolate module provides destroy_splines_* for its types. Every allocation in bspline_1d_init_uniform at line 43 has no corresponding deallocation path. Add a final(bspline_1d_finalize) procedure or a bspline_1d_destroy subroutine, and equivalent for 2D/3D.
  2. [blocker] src/interpolate/neo_bspline_interp.f90:54 — When dgbsv returns info /= 0 (singular matrix), the routine silently sets coeff = 0.0_dp and returns. The caller has no way to detect failure. Either signal through an intent(out) :: info argument, raise error stop, or document the failure mode. The same pattern occurs at lines 115, 121, 193, 203, 219 in the 2D/3D routines.
  3. [major] test/source/test_neo_bspline_2d.f90:77 — test_bspline_2d_interp_exact claims "direct interpolation (collocation)" with tolerance 1.0d-10, but evaluates cos(x1)*cos(2*x2) which is not a polynomial. The cubic spline space reproduces polynomials up to degree 3 exactly; cos is not one. The 3D test correctly uses a separable cubic polynomial. Either change the test function to a separable bivariate cubic (matching the 3D test pattern) or rename to test_bspline_2d_interp_approximate with a relaxed tolerance.
  4. [major] src/interpolate/neo_bspline_base.f90:44 — No input validation. If n_ctrl <= degree, h divides by zero or negative. If degree < 0 or n_ctrl < 1, the knot vector is invalid. Add minimal guards at the start of bspline_1d_init_uniform.
  5. [major] src/interpolate/neo_bspline_interp.f90:48 — The banded matrix column index col = span - p + j can be out of the 1:n range expected by dgbsv. When span = p+1 and j = 0, col = 1 (correct). When span = n and j = p, col = n+p-p = n (correct). However, the ab array is allocated as ab(ldab, n) with columns indexed 0:n-1 implicitly — Fortran arrays default to 1:n. The column col used at line 48 is 1-based from the B-spline control point indexing, but the ab allocation ab(ldab, n) creates indices 1:n. Verify that col never reaches 0 (it doesn't for uniform knots, but the code should be robust).
  6. [minor] .gitignore:69 — bench_neo_bspline_*.png is redundant; test/**/*.png on line 68 already covers this. Remove the redundant line.
  7. [minor] test/CMakeLists.txt:242 — bench_neo_bspline_vs_interpolate.x is compiled and linked but never registered with add_test. It will not run in CI. Either register it or remove it from this PR.
  8. [minor] src/interpolate/neo_bspline.f90:7 — Umbrella module uses use neo_bspline_base without only. This imports all public symbols from the submodule. Given the submodules declare private by default, this is safe but only: would make the re-export explicit and prevent accidental exposure if submodule visibility changes.
  9. [minor] src/interpolate/neo_bspline_base.f90:25 — find_span and basis_funs are public. These are implementation details of the Cox-de Boor algorithm. If they're needed only by interp and lsq submodules, make them module-procedure or keep them public only if downstream code needs direct basis access.
  10. [minor] test/source/test_neo_bspline_1d.f90:40 — Partition of unity test evaluates each basis function by constructing a full coefficient vector with one element set to 1.0. A more direct approach would sum basis_funs output directly, avoiding the N_CTRL evaluations per sample point. Not incorrect, just inefficient.

Verdict: Request changes — the memory leak (no finalize/destroy for allocatable knot vectors) and silent dgbsv failures are correctness issues that must be addressed before merge.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@slopqueue slopqueue[bot] slopqueue[bot] requested changes

Requested changes must be addressed to merge this pull request.

Assignees

No one assigned

Labels

Review effort 3/5

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

1 participant

@krystophny