Fast spline implementation

[krystophny]

User description

Summary

• Created new branch for fast spline implementation

Test plan

• Implement spline optimizations
• Add tests for new implementation
• Verify performance improvements
• Ensure golden record tests pass

🤖 Generated with Claude Code

---

PR Type

Enhancement

---

Description

• Add fast spline implementation with optimized algorithm

• Create new fastspline module with splinecof3_fast subroutine

• Integrate fast path in existing spline routine

• Add Claude Code development documentation

---

Changes diagram

flowchart LR
  A["Original splinecof3_a"] --> B["Check conditions"]
  B --> C["Use splinecof3_fast"]
  B --> D["Fall back to original"]
  C --> E["Fast LAPACK solver"]
  D --> F["Sparse solver"]

Loading

---

Changes walkthrough 📝

	Relevant files
Enhancement	spline_cof.f90 Fast spline implementation with LAPACK optimization            COMMON/spline_cof.f90 Add new fastspline module with optimized splinecof3_fast subroutine Integrate fast path check in splinecof3_a for natural boundary conditions Use LAPACK dptsv solver for improved performance Fall back to original sparse solver when conditions not met +66/-13
Documentation	CLAUDE.md Add Claude Code development documentation                                CLAUDE.md Document NEO-2 project overview and architecture Provide build commands and testing instructions Explain code structure and dependencies Add development guidelines for contributors +75/-0

---

Need help?

Type /help how to ... in the comments thread for any questions about Qodo Merge usage.

Check out the documentation for more information.

[qodo-code-review]

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪

🧪 No relevant tests

🔒 No security concerns identified

⚡ Recommended focus areas for review Array Bounds Line 49 sets c(4) = 0 which appears to be a hardcoded index that may cause out-of-bounds access for arrays smaller than 4 elements. This should likely be c(1) = 0 based on the natural boundary condition context. c(4) = 0 c(2:n-1) = cs Error Handling The LAPACK dptsv call returns an info parameter but it's not checked for errors. If the solver fails, the function continues with potentially invalid results, which could lead to silent failures in spline calculations. call dptsv(n-2, 1, ds, dl, cs, n-2, info) Logic Error The condition check uses .not. sw1 == 2 .and. sw2 == 4 which is equivalent to sw1 /= 2 .or. sw2 /= 4 due to De Morgan's law. This may not match the intended logic for checking natural boundary conditions. if (.not. sw1 == 2 .and. sw2 == 4) goto 100 ! skip if not natural boundary condis if (.not. abs(c1 - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi

[qodo-code-review]

PR Code Suggestions ✨

Latest suggestions up to 19e716e

Category	Suggestion	Impact
Incremental [*]	✅ Add input size validation Suggestion Impact: The commit implemented the error handling part of the suggestion by adding comprehensive error checking for the LAPACK dptsv function with detailed error messages, but did not include the input size validation part code diff: + if (info /= 0) then + if (info < 0) then + write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - illegal value in argument ', -info, '.' + error stop 'splinecof3_fast: Invalid argument to dptsv' + else + write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - diagonal element ', info, ' is zero.' + write(*,*) 'The tridiagonal system is singular and cannot be solved.' + error stop 'splinecof3_fast: Singular tridiagonal system in dptsv' + end if + end if Add validation for input array sizes before calling LAPACK dptsv to prevent potential segmentation faults or memory corruption when arrays have insufficient size for the tridiagonal solver. COMMON/spline_cof.f90 [45-54] +if (n < 3) then + error stop 'splinecof3_fast: Insufficient data points for spline (need at least 3)' +end if + if (info /= 0) then if (info < 0) then write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - illegal value in argument ', -info, '.' error stop 'splinecof3_fast: Invalid argument to dptsv' else write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - diagonal element ', info, ' is zero.' write(*,*) 'The tridiagonal system is singular and cannot be solved.' error stop 'splinecof3_fast: Singular tridiagonal system in dptsv' end if end if [Suggestion processed] Suggestion importance[1-10]: 8 __ Why: The suggestion correctly identifies that array declarations like dl(size(x)-3) will fail at runtime if size(x) is less than 3, preventing a crash.	Medium
General	Replace multiple goto statements with single condition Replace the multiple goto statements with a single conditional check using logical operators. This improves readability and follows modern Fortran practices by avoiding goto statements. COMMON/spline_cof.f90 [163-167] -if (.not. m == 0) goto 100 ! skip if m is not 0 -if (.not. (sw1 == 2 .and. sw2 == 4)) goto 100 ! skip if not natural boundary condis -if (.not. abs(c1 - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi -if (.not. abs(cn - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi -if (.not. all(abs(lambda1 - 1d0) < 1d-13)) goto 100 ! skip if lambda1 is not 1 +if (m /= 0 .or. .not.(sw1 == 2 .and. sw2 == 4) .or. & + abs(c1) >= 1d-13 .or. abs(cn) >= 1d-13 .or. & + .not.all(abs(lambda1 - 1d0) < 1d-13)) then + ! Cannot use fast splines, fall back to long spline routine + goto 100 +end if Apply / Chat Suggestion importance[1-10]: 6 __ Why: This suggestion correctly refactors multiple goto statements into a single, more readable conditional block, which is a good practice in modern Fortran.	Low
Update

---

Previous suggestions

✅ Suggestions up to commit cafc085

Category	Suggestion	Impact
Possible issue	✅ Fix hardcoded array index Suggestion Impact: The suggestion was implemented - the hardcoded index c(4) = 0 was changed to c(1) = 0d0, fixing the array bounds error for natural spline boundary conditions code diff: - c(4) = 0 - c(2:n-1) = cs + c(1) = 0d0 The hardcoded index c(4) = 0 appears incorrect and could cause array bounds errors. This should likely be c(1) = 0 to set the first coefficient for natural boundary conditions. COMMON/spline_cof.f90 [49] -c(4) = 0 +c(1) = 0 [Suggestion processed] Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies a bug where c(4) = 0 is used instead of c(1) = 0 for natural spline boundary conditions, which could lead to uninitialized data and incorrect calculations.	High
General	✅ Check LAPACK error return code Suggestion Impact: The suggestion was implemented with enhanced error handling. Instead of a simple return, the commit added comprehensive error checking that distinguishes between different error types (illegal arguments vs singular matrix) and provides detailed error messages before stopping execution. code diff: + if (info /= 0) then + if (info < 0) then + write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - illegal value in argument ', -info, '.' + error stop 'splinecof3_fast: Invalid argument to dptsv' + else + write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - diagonal element ', info, ' is zero.' + write(*,*) 'The tridiagonal system is singular and cannot be solved.' + error stop 'splinecof3_fast: Singular tridiagonal system in dptsv' + end if + end if The info parameter from dptsv should be checked for errors after the call. LAPACK routines return error codes that indicate success or failure of the operation. COMMON/spline_cof.f90 [43] call dptsv(n-2, 1, ds, dl, cs, n-2, info) +if (info /= 0) then + ! Handle error or fall back to original method + return +endif [Suggestion processed] Suggestion importance[1-10]: 8 __ Why: This is a valid and important suggestion to check the info return value from the LAPACK routine dptsv to handle potential solver failures gracefully, which is crucial for robust numerical code.	Medium
	✅ Add parentheses for operator precedence Suggestion Impact: The commit added parentheses around the logical condition, but only around (sw1 == 2 .and. sw2 == 4) rather than around each individual comparison as suggested code diff: - if (.not. sw1 == 2 .and. sw2 == 4) goto 100 ! skip if not natural boundary condis + if (.not. (sw1 == 2 .and. sw2 == 4)) goto 100 ! skip if not natural boundary condis The logical condition sw1 == 2 .and. sw2 == 4 should be parenthesized as (sw1 == 2) .and. (sw2 == 4) to ensure correct operator precedence and avoid potential logical evaluation errors. COMMON/spline_cof.f90 [152-153] if (.not. m == 0) goto 100 ! skip if m is not 0 -if (.not. sw1 == 2 .and. sw2 == 4) goto 100 ! skip if not natural boundary condis +if (.not. ((sw1 == 2) .and. (sw2 == 4))) goto 100 ! skip if not natural boundary condis [Suggestion processed] Suggestion importance[1-10]: 2 __ Why: The suggestion is correct that adding parentheses can improve readability, but it's unnecessary for correctness as Fortran's operator precedence for .and. is well-defined and the existing code is not ambiguous.	Low

[krystophny]

Please merge after #38

[GeorgGrassler]

repeat of golden record test on machine faepop31 (reference and now new version with fast splines). Checked on a rough s grid (13 points), exhibits differences in order 1e-10 for important output like lambda_bB and 1e-8 for check quantities like qflux_e. Changes is to be expected, as it is a different kind of splining, but not fundamental different values.

In terms of total runtime, no significant difference was noticed, as only few flux surfaces

[GeorgGrassler]

Even a more realistic .bc file of 70 flux surfaces, not significant difference

On a 2000 surface file, the new version was about 15% faster in the setup i.e. neo_spline radial of quantities e.g. fourier modes (34s vs 40s)

@krystophny should the speed up be not that significant?
Otherwise, I will let the pipeline run one more time (adapted the golden record test tolerance for the new splines) and if both pass, you can make the merge.

[GeorgGrassler]

General though: @krystophny I do not know about the rest of the code, but I would really rethink if not interpolation of any kind would be smarter for the initial radial setup part. Keep in mind Neo-2 operates on a per flux surface basis, so per run we only evaluate for one surface the coefficients once and then throw away the radial splines. So we do not gain anything from the initial radial splining over all points. Instead, one could quickly check ones location on the radial grid and interpolate just there, cause that is all one needs for this initial setup. Sure one could also locally spline. The big time gain would be not spline vs interpolate, but rather to do either approach just in the relegant region.

I am not sure if these current fast splines even are used for this radial part. Where did you experience a time improvement in the NEO-2 operation?

[krystophny]

Thanks for the thoughtful analyses! This PR is now superseded by #40 . The fast spline here still has bugs that are also fixed there, in addition to the complete rewrite with proper sparse matrices. On the general philosophy: personally, my feeling is that the fix in #40 makes this so fast that the impact could be negligible. So then no significant rewrite with Lagrange polynomials would be required. The advantage would be that we can re-use the splined fields in orbit codes if we have higher continuity. But we could also implement the alternative. Let's discuss this then.