Add low-level split-solve API for advanced users

[flaport]

Summary

Addresses #16 by exposing klu_analyze, klu_factor, and klu_solve separately for advanced users who need to optimize Newton-Raphson solvers by reusing symbolic analysis.

New Functions in klujax_cpp

True split-solve API (handle-based):

• klu_analyze(Bp, Bi) - Symbolic analysis on CSC sparsity pattern. Returns a handle.
• klu_factor_f64(symbolic_handle, Bp, Bi, Bx) - Numeric factorization (float64). Returns a handle.
• klu_factor_c128(symbolic_handle, Bp, Bi, Bx) - Numeric factorization (complex128). Returns a handle.
• klu_solve_f64(symbolic_handle, numeric_handle, b) - Solve using pre-computed factorizations (float64).
• klu_solve_c128(symbolic_handle, numeric_handle, b) - Solve using pre-computed factorizations (complex128).
• klu_free_symbolic(handle) - Free a symbolic factorization handle.
• klu_free_numeric(handle) - Free a numeric factorization handle.

FFI-based helpers:

• coo_to_csc - Convert COO format (Ai, Aj) to CSC format (Bp, Bi, Bk)
• solve_csc_f64 / solve_csc_c128 - Solve using pre-computed CSC structure

Usage Example (Newton-Raphson)

import numpy as np
import klujax_cpp

# COO to CSC conversion (one-time)
Bp, Bi, Bx = convert_coo_to_csc(...)  # your helper

# Symbolic analysis (one-time for fixed sparsity pattern)
sym_handle = klujax_cpp.klu_analyze(Bp, Bi)

# Newton-Raphson loop
for iteration in range(max_iterations):
    # Update Jacobian values (same sparsity, new values)
    Bx = compute_jacobian_values(...)
    
    # Numeric factorization (each iteration)
    num_handle = klujax_cpp.klu_factor_f64(sym_handle, Bp, Bi, Bx)
    
    # Solve
    dx = klujax_cpp.klu_solve_f64(sym_handle, num_handle, residual)
    
    # Free numeric handle (symbolic is reused!)
    klujax_cpp.klu_free_numeric(num_handle)
    
    x = x + dx
    if converged:
        break

# Cleanup
klujax_cpp.klu_free_symbolic(sym_handle)

The key benefit is that klu_analyze (the expensive symbolic analysis) only needs to be called once per sparsity pattern, while klu_factor can be called repeatedly with different values.

Note: These are low-level functions accessible via klujax_cpp module. The high-level klujax.solve API remains unchanged.

Test plan

• All 46 tests pass (including 5 new tests for split-solve API)
• New functions are exported and accessible via klujax_cpp
• Pre-commit hooks pass

Closes #16

🤖 Generated with Claude Code

[cdaunt]

Hi @flaport

Close, but not quite there yet :)

I will try and find some time to dig through in more detail but the quick feedback is that these do not appear to be jit compatible (and I did not see it in the test suite either)

All the tests use numpy and not jnp.numpy, is this intentional? For a non-linear+transitent solver the goal would be to pass jnp.arrays to a jitted function

`
File ~/code/circulus/circulus/solvers/strategies.py:366, in KLUSplitSolver._solve_impl(self, all_vals, residual)
    364     num_handle = klujax_cpp.klu_factor_c128(self.sym_handle, Bp, Bi, Bx)
    365 else:
--> 366     num_handle = klujax_cpp.klu_factor_f64(self.sym_handle, Bp, Bi, Bx)
    368 # 5. Solve
    369 if is_complex:

TypeError: klu_factor_f64(): incompatible function arguments. The following argument types are supported:
    1. (symbolic_handle: int, Bp: numpy.ndarray[numpy.int32], Bi: numpy.ndarray[numpy.int32], Bx: numpy.ndarray[numpy.float64]) -> int

Invoked with: 1, JitTracer<int64[11]>, JitTracer<int32[46]>, JitTracer<float64[46]>`

[flaport]

Closing in favor of #19