Use the truncated QR LAPACK routine geqp3rk as the backend for the partial qr factorization. fixes #60

[ajinkya-k]

The LAPACK routine geqpr3k available starting in v3.12.0 computes the truncated QR factorization based on one of three stopping conditions. This PR adds a backend for pqr that uses the LAPACK routine instead of the Julia native backend in the function geqp3_adap_main!. This backend is faster than using the current implementation.

There is still one unresolved issue that I need help with. When a max rank to truncate at is not provided, the R factor from the lapack backend differs from the current version and I am not sure why. Fixes #60

[ajinkya-k]

Regarding the discrepancy in the factor R, its because of a slightly different order of pivots. See more details in this comment in issue #60

[ajinkya-k]

relevant upstream issue in openBLAS OpenMathLib/OpenBLAS#5013

[ajinkya-k]

Hello! can someone help me with this PR? Since OpenBLAS was updated to v0.3.29 the LAPACK routines can be used for performing partial QR. Using the LAPACK routine is indeed faster. However, the tests fail for the Complex matrices. I dont see why that should be the case.

[MikaelSlevinsky]

Also, Julia 1.0 does not have the function require_one_based_indexing. This is a different error affecting the first half of the test environments.

[ajinkya-k]

@MikaelSlevinsky I still havent added compatibility. I am checking LAPACK version before invoking the LAPACK version. I only get test failures for complex 32 and complex 64. Here are some samples:

 none/ComplexF32
pqr: Test Failed at /Users/ajinkya/projects/julia_repos/testngdir/testing-pqr-lpk.jl:69
  Expression: norm(A - Matrix(F)) < approx_rtol * nrm
   Evaluated: 90.50967f0 < 0.0053947964f0

Stacktrace:
 [1] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:679 [inlined]
 [2] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:69 [inlined]
 [3] macro expansion
   @ ./timing.jl:581 [inlined]
 [4] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:44 [inlined]
 [5] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:1704 [inlined]
 [6] top-level scope
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:39
pqr: Test Failed at /Users/ajinkya/projects/julia_repos/testngdir/testing-pqr-lpk.jl:70
  Expression: norm(A' - Matrix(Fc)) < approx_rtol * nrm
   Evaluated: 90.50967f0 < 0.0053947964f0

Stacktrace:
 [1] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:679 [inlined]
 [2] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:70 [inlined]
 [3] macro expansion
   @ ./timing.jl:581 [inlined]
 [4] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:44 [inlined]
 [5] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:1704 [inlined]
 [6] top-level scope
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:39
pqr: Test Failed at /Users/ajinkya/projects/julia_repos/testngdir/testing-pqr-lpk.jl:74
  Expression: norm(y - F * xn) < approx_rtol * norm(y)
   Evaluated: 335.49106f0 < 0.019996826f0

Stacktrace:
 [1] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:679 [inlined]
 [2] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:74 [inlined]
 [3] macro expansion
   @ ./timing.jl:581 [inlined]
 [4] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:44 [inlined]
 [5] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:1704 [inlined]
 [6] top-level scope
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:39
pqr: Test Failed at /Users/ajinkya/projects/julia_repos/testngdir/testing-pqr-lpk.jl:75
  Expression: norm(y - F' * xm) < approx_rtol * norm(y)
   Evaluated: 500.39462f0 < 0.029825844f0

Stacktrace:
 [1] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:679 [inlined]
 [2] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:75 [inlined]
 [3] macro expansion
   @ ./timing.jl:581 [inlined]
 [4] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:44 [inlined]
 [5] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:1704 [inlined]
 [6] top-level scope
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:39

For complex 64:

  none/ComplexF64
pqr: Test Failed at /Users/ajinkya/projects/julia_repos/testngdir/testing-pqr-lpk.jl:69
  Expression: norm(A - Matrix(F)) < approx_rtol * nrm
   Evaluated: 90.50966799187782 < 1.0048591735576132e-11

Stacktrace:
 [1] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:679 [inlined]
 [2] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:69 [inlined]
 [3] macro expansion
   @ ./timing.jl:581 [inlined]
 [4] macro expansion
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:44 [inlined]
 [5] macro expansion
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Test/src/Test.jl:1704 [inlined]
 [6] top-level scope
   @ ~/projects/julia_repos/testngdir/testing-pqr-lpk.jl:39

[ajinkya-k]

The tests do pass for Float

[ajinkya-k]

I tracked the error to this Error: no BLAS/LAPACK library loaded for sgeqp3rk_64_(). This is because an old version of OpenBLAS is being used in the CI system. This was fixed recently and looks like a compatibility error. My code will only work with Julia 1.12+. The if statement LAPACK.version() < v"3.12.0" in pqr.jl should take care of it, but it is not. Maybe I should make it 3.12.1 and wait for it to be made available inside Julia.

[ajinkya-k]

I am testing this PR again with julia v1.12 and it now runs!
Moreover, all tests pass for both float64 and float32 types, but the approximation is off for complex types. I am trying to figure out if the geqp3rk routine is doing something drastically different from the pure julia implementation. This is an example of a discrepancy I get:

pqr: Test Failed at REPL[12]:39
  Expression: norm(y - transpose(F) * xm) < approx_rtol * norm(y)
   Evaluated: 565.17487f0 < 0.033687048f0

can anyone help me identify the issue?