Speed up SparseMatrixCSC(::AbstractQ) constructor and multiplication

[SobhanMP]

closes #194

[SobhanMP]

is there any reason to avoid force push in branches like this?

[codecov-commenter]

Codecov Report

Merging #196 (a556eeb) into main (134f67a) will increase coverage by 0.02%.
The diff coverage is 93.10%.

@@            Coverage Diff             @@
##             main     #196      +/-   ##
==========================================
+ Coverage   91.85%   91.87%   +0.02%     
==========================================
  Files          11       11              
  Lines        7118     7148      +30     
==========================================
+ Hits         6538     6567      +29     
- Misses        580      581       +1     

Impacted Files	Coverage Δ
src/solvers/spqr.jl	88.31% <50.00%> (-1.03%)	⬇️
src/sparsematrix.jl	95.33% <100.00%> (+0.04%)	⬆️
src/linalg.jl	84.58% <0.00%> (+0.08%)	⬆️
src/solvers/umfpack.jl	88.28% <0.00%> (+0.12%)	⬆️

📣 Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

[dkarrasch]

This doesn't look right. Q is a 4x4 matrix, but now you display some 4x2 matrix.

[SobhanMP]

after some looking up i think i found the problem. itertate of QRsparseQ is implemented by AbstractQ with lmul!. one lmul! per index to get the column and then take the right elemenent

by handling it on a per column basis, the run time of sparse(Q) for a random 100x100 matrix goes from 263.841 ms (20014 allocations: 763.41 MiB) to 5.943 ms (390189 allocations: 25.96 MiB)

[ViralBShah]

is there any reason to avoid force push in branches like this?

Force pushing in branches is fine.

[dkarrasch]

As for the original issue, I think we'll need to change

SparseArrays.jl/src/solvers/spqr.jl

Lines 306 to 307 in d367f11

	println(io, "Q factor:")
	show(io, mime, F.Q)

to something like

println(io, "Q factor: " * summary(F.Q))

or, alternatively, modify the display of QRSparseQ to

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, Q::QRSparseQ)
    summary(io, Q)
end

to only show the summary. Really, I think we shouldn't make it look like Qs are like arrays and shouldn't display any "elements".

[SobhanMP]

IMHO leaving _from_lmul!/_from_col as a function was ok.

[dkarrasch]

IMHO leaving _from_lmul!/_from_col as a function was ok.

What for? I don't think we should "deposit" unused internal functions just for fun. I also removed the sizehint option because that was performing the expensive work (multiplication) twice, so I don't think that would be worth it in any case.

[dkarrasch]

The current constructor does not apply to adjoints of AbstractQs, but would apply automatically in case JuliaLang/julia#46196 gets accepted. Also, the constructor turns out to be useful in multiplication. Previously, F.Q * F.R would fall back to _generic_matmul!! 😱 But now it makes use of the sparse constructor. 🎉

[SobhanMP]

IMHO leaving _from_lmul!/_from_col as a function was ok.

What for? I don't think we should "deposit" unused internal functions just for fun. I also removed the sizehint option because that was performing the expensive work (multiplication) twice, so I don't think that would be worth it in any case.

i agree with remove the sizehint but the function SparseMatrixCSC(::AbstractQ) isn't really a property of AbstractQs but anything that implement lmul!. it wouldn't have been an unused internal (at least not _from_lmul!) as it was called by SparseMatrixCSC.

[dkarrasch]

isn't really a property of AbstractQs but anything that implement lmul!

That's right, up to the restriction that additionally there is no direct/fast getindex and no iteration protocol. If we find more of such types, we can add them to the signature. Do you have something in mind already?

[SobhanMP]

apart from matrix free jacobians, no. also since this package has high inertia, maybe it's better to have an undocumented helper than adding it when needed?

[ViralBShah]

Off topic - We really do need to reduce the inertia on this package (and I'm doing my best to help speed things along). Some inertia is required since it is still part of Julia stdlibs and widely depended upon. What else can we do? Maybe best to open a new issue for discussion.

[dkarrasch]

Sounds good. Let's have a documented internal function then. 😉 And bikeshed the name! I'll make a proposal and then we see.

[SobhanMP]

looks good to me!