Support for Unitful sparse arrays? [feature request]

[JeffFessler]

The current code for lldl fails for Unitful sparse (or dense) arrays. MWE:

using SparseArrays: spdiagm
using LimitedLDLFactorizations: lldl
using Unitful: s
A = spdiagm(Float32.(1:3)*s^2) # unitful sparse array
lldl(A) # fails

ERROR: MethodError: no method matching lldl(::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Int64}, ::Type{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}})
Closest candidates are:
  lldl(::SparseArrays.SparseMatrixCSC{Tv, Ti}; kwargs...) where {Tv<:Number, Ti<:Integer} at ~/.julia/packages/LimitedLDLFactorizations/usbOW/src/LimitedLDLFactorizations.jl:493
  lldl(::SparseArrays.SparseMatrixCSC{Tv, Ti}, ::Type{Tf}; P, memory, droptol, α, α_increase_factor, check_tril) where {Tv<:Number, Ti<:Integer, Tf<:Real} at ~/.julia/packages/LimitedLDLFactorizations/usbOW/src/LimitedLDLFactorizations.jl:471

Some of the lldl code uses Number which is general enough to include unitful values, such as here: https://github.com/JuliaSmoothOptimizers/LimitedLDLFactorizations.jl/blob/51eec8dae729b9c14486ae3465edb700fb2130aa/src/LimitedLDLFactorizations.jl#L493

but then just a bit deeper into the code the types get restricted to Real, like here: https://github.com/JuliaSmoothOptimizers/LimitedLDLFactorizations.jl/blob/51eec8dae729b9c14486ae3465edb700fb2130aa/src/LimitedLDLFactorizations.jl#L480

I realize that Number is general enough to include complex number types and perhaps you don't want to support those? (Though that could be useful too, right?) My Unitful values are reinterpretable as real numbers, but currently not when stored in a sparse matrix: https://github.com/JuliaSparse/SparseArrays.jl/issues/289

I figure it's a long shot, but I thought I'd ask to see if there is any possibility of supporting more general Number types, especially Unitful values and maybe eventually complex values too? Since you have a pure Julia version, it seems like it would be possible!

[dpo]

@geoffroyleconte Would anything break if we changed Real to Number here?

[geoffroyleconte]

Probably not, but the factorization would not be accurate for complex matrices.

I opened https://github.com/JuliaSmoothOptimizers/LimitedLDLFactorizations.jl/pull/79. I left lldl with a Symmetric / Hermitian matrix with the Real type restriction, but I changed lldl(::SparseArrays.SparseMatrixCSC{Tv, Ti}; kwargs...) so that you can provide the matrix without the Symmetric wrapper.

However, with the example you provided I have this error

ERROR: DimensionError: 0.0f0 s^2 and 1.0f0 s^4 are not dimensionally compatible.
Stacktrace:
 [1] +(x::Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}  , y::Unitful.Quantity{Float32, 𝐓^4, U 
nitful.FreeUnits{(s^4,), 𝐓^4, nothing}} )
   @ Unitful C:\Users\Geoffroy Leconte\.julia\packages\Unitful\wRgqZ\src\quantities.jl:133
 [2] macro expansion
   @ C:\Users\Geoffroy Leconte\.julia\dev\LimitedLDLFactorizations\src\LimitedLDLFactorizations.jl:303 [inlined]
 [3] macro expansion
   @ .\simdloop.jl:77 [inlined]
 [4] lldl_factorize!(S::LimitedLDLFactorizations.LimitedLDLFactorization{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnit 
s{(s^2,), 𝐓^2, nothing}}, Int64, Vector{Int64}, Vector{Int64}} , T::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}, Int64}  ; droptol::Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^ 
2,), 𝐓^2, nothing}} )
   @ LimitedLDLFactorizations C:\Users\Geoffroy Leconte\.julia\dev\LimitedLDLFactorizations\src\LimitedLDLFactorizations.jl:301
 [5] lldl(A::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}, Int6  
4}, ::Type{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}}  ; P::Vector{Int64}, memory::Int64, droptol::Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}  , α::Int64, α_increase_factor::Int64, check_tril::Bool)
   @ LimitedLDLFactorizations C:\Users\Geoffroy Leconte\.julia\dev\LimitedLDLFactorizations\src\LimitedLDLFactorizations.jl:490
 [6] lldl(A::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}, Int6  
4}, ::Type{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}}  )
   @ LimitedLDLFactorizations C:\Users\Geoffroy Leconte\.julia\dev\LimitedLDLFactorizations\src\LimitedLDLFactorizations.jl:471
 [7] lldl(A::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}, Int6  
4}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ LimitedLDLFactorizations C:\Users\Geoffroy Leconte\.julia\dev\LimitedLDLFactorizations\src\LimitedLDLFactorizations.jl:493
 [8] lldl(A::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float32, 𝐓^2, Unitful.FreeUnits{(s^2,), 𝐓^2, nothing}}, Int6  
4})
   @ LimitedLDLFactorizations C:\Users\Geoffroy Leconte\.julia\dev\LimitedLDLFactorizations\src\LimitedLDLFactorizations.jl:493
 [9] top-level scope
   @ REPL[17]:1

[JeffFessler]

Thanks for this start!

The MWE in #78 gets further with the code in this PR, but as you noted it fails on this line: https://github.com/geoffroyleconte/LimitedLDLFactorizations.jl/blob/2d747e8d5065e8855a250f53c821e6ec33e61f27/src/LimitedLDLFactorizations.jl#L303 in this block of code:

@inbounds @simd for col = 1:n
    dpcol = adiag[P[col]]
    wa1[col] += dpcol * dpcol
    wa1[col] = sqrt(wa1[col])
    wa1[col] > 0 && (s[col] = 1 / sqrt(wa1[col]))
  end

That code sequence is perplexing for data with units because the same location wa1[col] would be asked to store something with the square of the units of adiag in the += line and then back to the original units in the sqrt line, which is impossible.

If A has units s^2, like in my MWE, then in a LL' factorization the L must have units s. Here it is LDL' so I guess D could have the same units as A and L could be unitless, or D could be unitless and L would have the sqrt units. Based on the code block above, I think probably having D be unitless is easier? So i think the LLDL type might need to be generalized to have something like Ta for the eltype of A and Tl for the eltype of L where Tl = eltype(sqrt(oneunit(Ta))). I've gone through that type of unit analysis for repos that I manage for algorithms that I understand. Here I can help but don't know enough to provide a PR myself.

I modified my MWE as follows to use the Tf positional argument of lldl like this:

using SparseArrays: spdiagm
using LimitedLDLFactorizations: lldl
using Unitful: s
A = spdiagm(Float32.(1:3)*s^2) # unitful sparse array
Ta = eltype(A)
Tl = eltype(sqrt(oneunit(Ta)))
lldl(A, Tl) # fails

That version fails on this line: https://github.com/geoffroyleconte/LimitedLDLFactorizations.jl/blob/2d747e8d5065e8855a250f53c821e6ec33e61f27/src/LimitedLDLFactorizations.jl#L132

I suspect that small modifications of the types in the struct might do it...