UnitfulFactors.jl

using SparseArrays, LinearAlgebra, SuiteSparse
using Unitful

"""
    elunit(x)
    
Returns the unit of the elements of `x`.
(Returns a `Unitful.FreeUnits{(),NoDims,nothing}` object `x`'s unit is heterogeneous.)
"""
elunit(::SparseMatrixCSC{U,Int64}) where {U<:Quantity} = unit(U)
elunit(::AbstractVecOrMat{U}) where {U<:Quantity} = unit(U)

Base.:*(M::SparseMatrixCSC{<:Quantity,Int64}, x::VecOrMat{<:Quantity}) = (elunit(x) * elunit(M)) * (ustrip.(M) * ustrip.(x)) 

"""
    UnitfulFactors{T,U}
    
Type for LinearAlgebra to work with Unitful.
"""
struct UnitfulFactors{T,U}
  factors::T
  unit::U
end

LinearAlgebra.factorize(M::SparseMatrixCSC{<:Quantity,Int64}) = UnitfulFactors(factorize(ustrip.(M)), elunit(M))
LinearAlgebra.factorize(M::Array{<:Quantity,2}) = UnitfulFactors(factorize(ustrip.(M)), elunit(M))
LinearAlgebra.:\(M::SparseMatrixCSC{<:Quantity,Int64}, x::VecOrMat{<:Quantity}) = (ustrip.(M) \ ustrip.(x)) * (elunit(x) / elunit(M))
LinearAlgebra.:\(M::Array{<:Quantity,2}, x::VecOrMat{<:Quantity}) = (ustrip.(M) \ ustrip.(x)) * (elunit(x) / elunit(M))
LinearAlgebra.:\(F::UnitfulFactors, x::VecOrMat{<:Quantity}) = (F.factors \ ustrip.(x)) * (elunit(x) / F.unit)

# Tests

T = (sprand(10,10,0.5) + I) * u"1/s" ;
x = rand(10) * u"mmol/m^3" ;

T * x
T \ (T * x)
F = factorize(T)
F \ (T * x)

M = Matrix(ustrip.(T)) * elunit(T) ;
M * x
M \ (M * x)
F = factorize(M)
F \ (M * x)