Levenberg-Marquardt with Sparse J: J isn't square?!

[luke-kiernan]

Describe the bug 🐞

I'm trying to run the LevenbergMarquardt solver on the system of equations arising from a power flow problem. After setting everything up, NonlinearSolve.jl is saying that the Jacobian matrix isn't square. All matrices I've provided are indeed square.

Minimal Reproducible Example 👇

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)

Error & Stacktrace ⚠️

ERROR: LoadError: ArgumentError: KLU only accepts square matrices.
Stacktrace:
  [1] LinearSolveSparseArraysExt.KLU.KLUFactorization(A::SparseMatrixCSC{Float64, Int32})
    @ LinearSolveSparseArraysExt.KLU ~/.julia/packages/LinearSolve/O6GSm/src/KLU/klu.jl:177
  [2] init_cacheval
    @ ~/.julia/packages/LinearSolve/O6GSm/ext/LinearSolveSparseArraysExt.jl:145 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:308 [inlined]
  [4] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::SparseMatrixCSC{Float64, Int32}, b::Vector{Float64}, u::Vector{Float64}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool, assump::OperatorAssumptions{Bool})
    @ LinearSolve ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:291
  [5] #init#3
    @ ~/.julia/packages/LinearSolve/O6GSm/src/common.jl:239 [inlined]
  [6] init
    @ ~/.julia/packages/LinearSolve/O6GSm/src/common.jl:140 [inlined]
  [7] #init#95
    @ ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:270 [inlined]
  [8] init
    @ ~/.julia/packages/LinearSolve/O6GSm/src/default.jl:266 [inlined]
  [9] #construct_linear_solver#56
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/linear_solve.jl:77 [inlined]
 [10] construct_linear_solver
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/linear_solve.jl:61 [inlined]
 [11] init(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, alg::DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, J::SparseMatrixCSC{Float64, Int32}, fu::Vector{Float64}, u::Vector{Float64}; stats::SciMLBase.NLStats, pre_inverted::Val{false}, linsolve_kwargs::@NamedTuple{abstol::Float64, reltol::Float64}, abstol::Float64, reltol::Float64, timer::Nothing, alias_J::Bool, shared::Val{2}, kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/yZeYz/src/descent/damped_newton.jl:139
 [12] init
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/descent/damped_newton.jl:44 [inlined]
 [13] init(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, alg::GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, J::SparseMatrixCSC{Float64, Int32}, fu::Vector{Float64}, u::Vector{Float64}; shared::Val{1}, pre_inverted::Val{false}, linsolve_kwargs::@NamedTuple{abstol::Float64, reltol::Float64}, abstol::Float64, reltol::Float64, internalnorm::typeof(NonlinearSolveBase.L2_NORM), kwargs::@Kwargs{stats::SciMLBase.NLStats, timer::Nothing})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/yZeYz/src/descent/geodesic_acceleration.jl:83
 [14] __init(::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion{Float64}, GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, Nothing, Nothing, Nothing, Val{true}}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Nothing, reltol::Nothing, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
    @ NonlinearSolveFirstOrder ~/.julia/packages/NonlinearSolveFirstOrder/UpLMt/src/solve.jl:167
 [15] __init
    @ ~/.julia/packages/NonlinearSolveFirstOrder/UpLMt/src/solve.jl:121 [inlined]
 [16] #__solve#146
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/solve.jl:5 [inlined]
 [17] __solve
    @ ~/.julia/packages/NonlinearSolveBase/yZeYz/src/solve.jl:1 [inlined]
 [18] #solve_call#36
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:667 [inlined]
 [19] solve_call
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:624 [inlined]
 [20] solve_up(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion{Float64}, GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, Nothing, Nothing, Nothing, Val{true}}; kwargs::@Kwargs{alias_u0::Bool})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1183
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1177 [inlined]
 [22] #solve#44
    @ ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1171 [inlined]
 [23] solve(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, SciMLBase.AutoSpecialize, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, typeof(jac!), Nothing, Nothing, SparseMatrixCSC{Float64, Int32}, SparseMatrixCSC{Float64, Int32}, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, args::GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion{Float64}, GeodesicAcceleration{DampedNewtonDescent{Nothing, Float64, NonlinearSolveFirstOrder.LevenbergMarquardtDampingFunction{Float64, Float64, Float64}}, Float64, Float64}, Nothing, Nothing, Nothing, Val{true}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Qw6nY/src/solve.jl:1136
 [24] top-level scope
    @ ~/Documents/julia/nonlinear-solve-mwe/mwe.jl:37

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Status `~/Documents/julia/nonlinear-solve-mwe/Project.toml`
  [8913a72c] NonlinearSolve v4.8.0
  [94fada2c] PowerFlows v0.9.0
  [f00506e0] PowerSystemCaseBuilder v1.3.11
  [bcd98974] PowerSystems v4.6.2
  [2f01184e] SparseArrays v1.11.0

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

Status `~/Documents/julia/nonlinear-solve-mwe/Manifest.toml`
  [47edcb42] ADTypes v1.14.0
  [1520ce14] AbstractTrees v0.4.5
  [7d9f7c33] Accessors v0.1.42
  [79e6a3ab] Adapt v4.3.0
  [13e28ba4] AppleAccelerate v0.4.0
⌃ [4fba245c] ArrayInterface v7.18.0
  [4c555306] ArrayLayouts v1.11.1
  [62783981] BitTwiddlingConvenienceFunctions v0.1.6
  [a74b3585] Blosc v0.7.3
  [70df07ce] BracketingNonlinearSolve v1.2.0
  [2a0fbf3d] CPUSummary v0.2.6
  [336ed68f] CSV v0.10.15
  [d360d2e6] ChainRulesCore v1.25.1
  [fb6a15b2] CloseOpenIntervals v0.1.13
  [944b1d66] CodecZlib v0.7.8
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.1
  [f70d9fcc] CommonWorldInvalidations v1.0.0
  [34da2185] Compat v4.16.0
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.8
  [adafc99b] CpuId v0.3.1
  [a8cc5b0e] Crayons v4.1.1
  [a10d1c49] DBInterface v2.6.1
  [9a962f9c] DataAPI v1.16.0
  [a93c6f00] DataFrames v1.7.0
  [864edb3b] DataStructures v0.18.22
  [e2d170a0] DataValueInterfaces v1.0.0
  [8bb1440f] DelimitedFiles v1.9.1
⌃ [2b5f629d] DiffEqBase v6.170.1
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.6.52
  [ffbed154] DocStringExtensions v0.9.4
  [4e289a0a] EnumX v1.0.5
  [f151be2c] EnzymeCore v0.8.8
  [e2ba6199] ExprTools v0.1.10
  [55351af7] ExproniconLite v0.10.14
  [7034ab61] FastBroadcast v0.3.5
  [9aa1b823] FastClosures v0.3.2
  [a4df4552] FastPower v1.1.2
  [48062228] FilePathsBase v0.9.24
  [1a297f60] FillArrays v1.13.0
  [6a86dc24] FiniteDiff v2.27.0
  [f6369f11] ForwardDiff v1.0.1
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [46192b85] GPUArraysCore v0.2.0
  [c8ec2601] H5Zblosc v0.1.2
  [f67ccb44] HDF5 v0.17.2
  [615f187c] IfElse v0.1.1
  [2cd47ed4] InfrastructureSystems v2.6.0
  [842dd82b] InlineStrings v1.4.3
  [3587e190] InverseFunctions v0.1.17
  [41ab1584] InvertedIndices v1.3.1
  [92d709cd] IrrationalConstants v0.2.4
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.7.0
  [0f8b85d8] JSON3 v1.14.2
  [ae98c720] Jieko v0.2.1
  [ef3ab10e] KLU v0.6.0
  [ba0b0d4f] Krylov v0.10.1
  [b964fa9f] LaTeXStrings v1.4.0
  [10f19ff3] LayoutPointers v0.1.17
  [5078a376] LazyArrays v2.6.1
  [1d6d02ad] LeftChildRightSiblingTrees v0.2.0
  [87fe0de2] LineSearch v0.1.4
⌃ [7ed4a6bd] LinearSolve v3.9.0
  [2ab3a3ac] LogExpFunctions v0.3.29
⌅ [33e6dc65] MKL v0.7.0
  [3da0fdf6] MPIPreferences v0.1.11
  [1914dd2f] MacroTools v0.5.16
  [d125e4d3] ManualMemory v0.1.8
  [bb5d69b7] MaybeInplace v0.1.4
  [e1d29d7a] Missings v1.2.0
  [2e0e35c7] Moshi v0.3.5
  [46d2c3a1] MuladdMacro v0.2.4
  [ffc61752] Mustache v1.0.20
  [77ba4419] NaNMath v1.1.3
  [8913a72c] NonlinearSolve v4.8.0
  [be0214bd] NonlinearSolveBase v1.6.0
  [5959db7a] NonlinearSolveFirstOrder v1.4.0
  [9a2c21bd] NonlinearSolveQuasiNewton v1.3.0
  [26075421] NonlinearSolveSpectralMethods v1.2.0
  [bac558e1] OrderedCollections v1.8.0
  [d96e819e] Parameters v0.12.3
⌅ [46dd5b70] Pardiso v0.5.7
  [69de0a69] Parsers v2.8.3
  [f517fe37] Polyester v0.7.16
  [1d0040c9] PolyesterWeave v0.2.2
  [2dfb63ee] PooledArrays v1.4.3
  [dd99e9e3] PowerFlowData v1.5.0
  [94fada2c] PowerFlows v0.9.0
⌅ [bed98974] PowerNetworkMatrices v0.12.1
  [f00506e0] PowerSystemCaseBuilder v1.3.11
  [bcd98974] PowerSystems v4.6.2
⌅ [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [08abe8d2] PrettyTables v2.4.0
  [33c8b6b6] ProgressLogging v0.1.4
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.33.0
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.14
  [94e857df] SIMDTypes v0.1.0
  [0aa819cd] SQLite v1.6.1
  [0bca4576] SciMLBase v2.86.2
  [19f34311] SciMLJacobianOperators v0.1.3
⌃ [c0aeaf25] SciMLOperators v0.3.13
  [53ae85a6] SciMLStructures v1.7.0
  [91c51154] SentinelArrays v1.4.8
  [efcf1570] Setfield v1.1.2
⌃ [727e6d20] SimpleNonlinearSolve v2.3.0
  [a2af1166] SortingAlgorithms v1.2.1
  [0a514795] SparseMatrixColorings v0.4.19
  [276daf66] SpecialFunctions v2.5.1
  [aedffcd0] Static v1.2.0
  [0d7ed370] StaticArrayInterface v1.8.0
  [1e83bf80] StaticArraysCore v1.4.3
  [10745b16] Statistics v1.11.1
  [7792a7ef] StrideArraysCore v0.5.7
  [69024149] StringEncodings v0.3.7
  [892a3eda] StringManipulation v0.4.1
  [856f2bd8] StructTypes v1.11.0
  [2efcf032] SymbolicIndexingInterface v0.3.40
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [5d786b92] TerminalLoggers v0.1.7
  [8290d209] ThreadingUtilities v0.5.3
  [9e3dc215] TimeSeries v0.24.2
  [a759f4b9] TimerOutputs v0.5.28
  [3bb67fe8] TranscodingStreams v0.11.3
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [ea10d353] WeakRefStrings v1.4.2
  [76eceee3] WorkerUtilities v1.6.1
  [ddb6d928] YAML v0.4.13
  [0b7ba130] Blosc_jll v1.21.7+0
  [0234f1f7] HDF5_jll v1.14.6+0
  [e33a78d0] Hwloc_jll v2.12.0+0
⌅ [1d5cc7b8] IntelOpenMP_jll v2024.2.1+0
  [94ce4f54] Libiconv_jll v1.18.0+0
  [5ced341a] Lz4_jll v1.10.1+0
⌅ [856f044c] MKL_jll v2024.2.0+0
  [7cb0a576] MPICH_jll v4.3.0+1
  [f1f71cc9] MPItrampoline_jll v5.5.3+0
  [9237b28f] MicrosoftMPI_jll v10.1.4+3
  [fe0851c0] OpenMPI_jll v5.0.7+2
  [458c3c95] OpenSSL_jll v3.5.0+0
  [efe28fd5] OpenSpecFun_jll v0.5.6+0
  [76ed43ae] SQLite_jll v3.48.0+0
  [3161d3a3] Zstd_jll v1.5.7+1
  [477f73a3] libaec_jll v1.1.3+0
  [1317d2d5] oneTBB_jll v2022.0.0+0
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [05823500] OpenLibm_jll v0.8.5+0
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`

• Output of versioninfo()

Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 12 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

Additional context

There's a chance that I'm just doing something wrong here, that this isn't actually a bug. But the error message is bizarre enough that I decided to pick "Bug" over "Question" when opening this issue. I also created a discourse post about this same issue.

[avik-pal]

LM needs to solve a linear system with a size 2N x N matrix. The bottom N x N part is a diagonal term. I am not sure why we are dispatching to KLU, so this might be a LinearSolve issue?

[ChrisRackauckas]

https://github.com/SciML/LinearSolve.jl/blob/main/ext/LinearSolveSparseArraysExt.jl#L248

Can someone add an @show into there and see if assump.issq is true?

[luke-kiernan]

Can someone add an @show into there and see if assump.issq is true?

I did so, and I'm getting that it's false. I also added a print statement to line 145 of the same file: indeed, the size of A is (20,10).

Edit: working my way up the stacktrace, the errant matrix is created in stack frame [12] or higher. i.e. if I add a print statement around line 138 of damped_newton.jl inside NonlinearSolveBase, the matrix A has wrong dimensions there too.

[ChrisRackauckas]

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")

ERROR: KeyError: key :PowerSystems not found
Stacktrace:
  [1] root_module(where::Module, name::Symbol)
    @ Base ./loading.jl:1951
  [2] get_module(module_name::String)
    @ InfrastructureSystems ~/.julia/packages/InfrastructureSystems/yUXE4/src/utils/utils.jl:382
  [3] get_type_from_serialization_metadata(metadata::Dict{String, Any})
    @ InfrastructureSystems ~/.julia/packages/InfrastructureSystems/yUXE4/src/serialization.jl:132
  [4] get_type_from_serialization_data(data::Dict{String, Any})
    @ InfrastructureSystems ~/.julia/packages/InfrastructureSystems/yUXE4/src/serialization.jl:128
  [5] deserialize_components!(sys::PowerSystems.System, raw::Dict{String, Any})
    @ PowerSystems ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:1464
  [6] deserialize(::Type{…}, filename::String; time_series_read_only::Bool, time_series_directory::Nothing, kwargs::@Kwargs{…})
    @ PowerSystems ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:1441
  [7] deserialize
    @ ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:1389 [inlined]
  [8] PowerSystems.System(file_path::String; assign_new_uuids::Bool, kwargs::@Kwargs{})
    @ PowerSystems ~/.julia/packages/PowerSystems/9ugJZ/src/base.jl:204
  [9] _build_system(category::Type{…}, name::String, sys_descriptor::SystemDescriptor, print_stat::Bool; kwargs::@Kwargs{})
    @ PowerSystemCaseBuilder ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:112
 [10] _build_system
    @ ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:75 [inlined]
 [11] #build_system#236
    @ ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:26 [inlined]
 [12] build_system
    @ ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:18 [inlined]
 [13] build_system(category::Type{MatpowerTestSystems}, name::String)
    @ PowerSystemCaseBuilder ~/.julia/packages/PowerSystemCaseBuilder/7w4lG/src/build_system.jl:18
 [14] top-level scope
    @ ~/Desktop/test2.jl:8
Some type information was truncated. Use `show(err)` to see complete types.

[ChrisRackauckas]

And I couldn't recreate it with a simple MWE:

using NonlinearSolve

function f!(du, u, p)
    du[1] = 2u[1] - 2
    du[2] = u[1] - 4u[2]
end

jac = SparseMatrixCSC([
    1 0
    1 1
])

u0 = [0.0, 0.0]
prob = NonlinearProblem(
    NonlinearFunction(f!, jac_prototype = jac), u0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)

Please try and make an MWE that doesn't require the extra packages. The core issue shouldn't require anything specific about f, it's just how you setup the problem, so it would be good to just choose a simple f and recreate the rest.

[luke-kiernan]

ERROR: KeyError: key :PowerSystems not found

Try adding PowerSystems as well.

And I couldn't recreate it with a simple MWE:

Strangely, if I change the indexing type to Int32, then I get the same error:

using NonlinearSolve
using SparseArrays

function f!(du, u, p)
    du[1] = 2u[1] - 2
    du[2] = u[1] - 4u[2]
end

jac = SparseMatrixCSC{Float64, Int32}([
    1.0 0.0
    1.0 1.0
])

u0 = [0.0, 0.0]
prob = NonlinearProblem(
    NonlinearFunction(f!, jac_prototype = jac), u0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)

[ChrisRackauckas]

Thanks, that pinpointed where the issue was. ArrayInterface and LinearSolve need to handle Int32 cases, so it'll need to wait for the LinearSolve.jl bump but this is fixed locally and just needs one more merge.