# NonlinearSolve.jl with Sparse Jacobian: J isn't square?!

**URL:** https://discourse.julialang.org/t/nonlinearsolve-jl-with-sparse-jacobian-j-isnt-square/128768
**Category:** General Usage
**Tags:** nonlinearsolve
**Created:** [May 6, 2025, 8:22pm UTC](https://discourse.julialang.org/t/nonlinearsolve-jl-with-sparse-jacobian-j-isnt-square/128768 "2025-05-06T20:22:15Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![luke-kiernan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luke-kiernan/32/214309_2.png) [@luke-kiernan](https://discourse.julialang.org/u/luke-kiernan)
#### Post date: [May 6, 2025, 8:22pm UTC](https://discourse.julialang.org/t/nonlinearsolve-jl-with-sparse-jacobian-j-isnt-square/128768/1 "2025-05-06T20:22:15Z")

</div>

I’m trying to use `NonlinearSolve.jl`’s `LevenbergMarquardt` solver on a powerflow. The `solve` call yields an error: `LoadError: ArgumentError: KLU only accepts square matrices`. That seems nonsensical: `J0` is square. I expect things to have dimension equal to twice the number of buses in the system: in this case, that’ll be 10. The print statements I’ve added bear that out. What’s the issue?

Possibly relevant: [sparsity detection](https://docs.sciml.ai/NonlinearSolve/stable/basics/sparsity_detection/)

```julia
using NonlinearSolve
using SciMLBase
using SparseArrays
using PowerFlows, PowerSystemCaseBuilder
const PF = PowerFlows
const PSB = PowerSystemCaseBuilder
# small system for demo purposes: actual target system has 10k buses. 
sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_case5_sys")
pf = PF.ACPowerFlow()
data = PF.PowerFlowData(pf, sys)
residual = PF.ACPowerFlowResidual(data, 1)
J = PF.ACPowerFlowJacobian(data, 1) # this initializes the sparse structure.
J0 = deepcopy(J.Jv) # take a copy of the sparse structure.
x0 = PF.calculate_x0(data, 1)
@assert size(J0, 1) == size(J0, 2) # yes, dimensions match. 
@assert size(J0, 1) == size(x0, 1)

function f!(du::Vector{Float64}, u::Vector{Float64}, ::Any)
    println(size(du))
    println(size(u))
    residual(du, u, 1) # in-place compute f at u
end

function jac!(A::SparseMatrixCSC{Float64, Int32}, u::Vector{Float64}, ::Any)
    println(size(u))
    println(size(A))
    PF.update_data!(data, u, 1)
    J(A, 1) # in-place compute J at the point stored in data (set to u above)
end

func = NonlinearFunction{true, SciMLBase.AutoSpecialize}(
    f!;
    jac = jac!,
    jac_prototype = J0,
)

prob = NonlinearProblem(func, x0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)

```