Make QRPivotedWs compatible with ormqr!

[cohensbw]

I was thrilled when I discovered your package, exactly what I was looking for!
Great docs and super easy to use, I definitely plan on using it in some of my future projects!

For one of the project I am currently working on I need to call geqp3! and then ormqr! (or orgqr! instead) afterwards to construct construct the Q matrix.
How challenging would it be to make it so that the QRPivotedWs workspace is compatible with ormqr! (or orgqr!)?

One idea would also be make something like abstract type QRWorkspace <: Workspace end, and then create a version of ormqr! (or orgqr!) that take any inherited type from QRWorkspace as an argument?

Just wanted to say again, I really like your package and how well you have documented it, so thank you for developing it!

[louisponet]

Hi,

Fantastic that the package is of use to you!

If what I found here is the correct way, then I think that should be very straightforward and I'll implement it shortly.

I'm not that well versed in using these methods myself, would the Compact QRWYWs also work similarly? Or would it use a different lapack function?

Cheers,
Louis

[cohensbw]

Thank you for responding and sorry about taking so long to respond, I have been traveling.
Yes, I believe the QRWYWs would also work similarly!
And the orgqr! method should also behave essentially the same for all the QR decomposition variants as well.

[louisponet]

Hi again,

I think #24 should have taken care of everything here. I'm not sure that QRWYWs would have worked too since it uses a matrix rather than a vector for the representation of Q. If not, feel free to open a PR with the right implementation.

Thanks again for the interest and input, it's very valuable to us!

[cohensbw]

Thank you so much for doing that, and I am sorry it took me so long to respond! I will test it out, but looking at what you did I think this is exactly what I need!

[cohensbw]

Your implementation works beautifully for real numbers.
Unfortunately, an error message is thrown when a complex number is passed.
In the case that B is a complex number, if I run the command

QRws = QRPivotedWs(B)

I get the following error message:

ReadOnlyMemoryError()

Stacktrace:
 [1] resize!(ws::QRPivotedWs{ComplexF64}, A::Matrix{ComplexF64}; work::Bool)
   @ FastLapackInterface ~/.julia/packages/FastLapackInterface/46eAD/src/qr.jl:278
 [2] resize!
   @ ~/.julia/packages/FastLapackInterface/46eAD/src/qr.jl:269 [inlined]
 [3] QRPivotedWs(A::Matrix{ComplexF64})
   @ FastLapackInterface ~/.julia/packages/FastLapackInterface/46eAD/src/qr.jl:287
 [4] eval
   @ ./boot.jl:368 [inlined]
 [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1428

I believe an additional temporary array is required when calling geqp3! for complex numbers.

[louisponet]

Interesting, I will investigate.

[cohensbw]

Thanks for fixing that!
I wanted to share with you the StableLinearAlgebra.jl package you assistance helped me develop (just the dev docs work right now, not the stable ones, as I haven't registered the package yet):
https://github.com/cohensbw/StableLinearAlgebra.jl
The algorithms exported by the package are essential in writing something called a determinant quantum Monte Carlo code of simulating interacting electrons.

[MichelJuillard]

Thanks @cohensbw I'm happy that our package is useful to you! I'm curious why do you need numerically stable algorithms in your field: very big matrices? A mix of very small and very large quantities?

[cohensbw]

Your second guess it essentially correct. In determinant quantum Monte Carlo (DQMC) simulations we need to multiply long chains of matrices together, and information about the small numbers is lost, and that information is important to retain, as we need to invert the sums of these long matrix products.

There is actually another Julia package very (very) similar to mine that exports almost the same functionality:
https://github.com/carstenbauer/StableDQMC.jl
This other package has a very nice SciPost paper associated with it as well:
https://scipost.org/10.21468/SciPostPhysCore.2.2.011

However, for various reasons (including eliminating dynamic memory allocations) I decided to implement my own version.

[MichelJuillard]

Excellent! Thanks for the background information.

[louisponet]

That is fantastic @cohensbw, I'm a physicist myself working on stuff not too far removed from your field (i'm doing magnetism, multiferroics, etc. from first principles).