Add wrapper for ?geqrfp

[mdhaber]

See NumFOCUS Small Development Grant Proposal 2019 Round 3 proposal for background.

Add wrapper for ?geqrfp.

Existing wrapper for ?geqrf provides some guidance.

Suggested signature:

! rq_a,tau,work,info = geqrfp(a,lwork=3*n,overwrite_a=0)

We will also want to create a routine to compute the optimal lwork value:

! work,info = geqrfp_lwork(n)

I'm a bit confused about how overwrite_a works - is that optional argument added automatically by F2PY? It is an optional argument for ?geqrf but there's no "code" for it in the wrapper.
Update: See @ilayn's response below.

If the user provides lwork=-1, it seems that ?geqrfp is supposed to determine the optimal lwork rather than actually solving the problem. I suppose the idea is that if you don't know what to choose for lwork, you call ?geqrfp once with lwork=-1 just to get the optimal value of lwork, then call it again with that value of lwork to actually solve the problem. Is that right? Should the wrapper go through that process if no lwork is provided by the user? The wrapper for ?geqrf doesn't use this approach, it seems; it just defaults to 3n. Where does 3n come from? In the NAG manual, it looks like 64*n is used.
Update: See @ilayn's response below.

As for tests, the NAG manual provides an example for ?geqrf, so these could be adapted for ?geqrfp. Alternatively, we could follow the approach used by scipy.linalg.qr to convert the output into matrices Q and R and write tests based on those, e.g.

• Create random matrix A
• Decompose using wrapper
• Convert to Q and R using approach in scipy.linalg.qr
• Verify that R is upper triangular and has all non-negative diagonal
• Verify that Q is orthogonal
• Verify that A = Q @ R

etc.
Probably should do this for an example where regular QR gives negative entries on diagonal of R. We should also check for empty matrix and both tall and wide matrices as input.

[ilayn]

I'm a bit confused about how overwrite_a works - is that optional argument added automatically by F2PY? It is an optional argument for ?geqrf but there's no "code" for it in the wrapper.

When variable declaration has the option copy in the intent it is generated automatically by f2py

If the user provides lwork=-1, it seems that ?geqrfp is supposed to determine the optimal lwork rather than actually solving the problem. I suppose the idea is that if you don't know what to choose for lwork, you call ?geqrfp once with lwork=-1 just to get the optimal value of lwork, then call it again with that value of lwork to actually solve the problem. Is that right? Should the wrapper go through that process if no lwork is provided by the user?

Calling the same routine twice (one with lwork=-1) is a decades old trick. The main reason is that a variable, NB, changes from hardware to hardware for optimal size computation.

The user could call with -1, but what we do is to wrap another routine <lapack_name>_lwork(...) (see examples in the wrapper files) so that you call that and receive the lwork number. Then you call the actual routine with that lwork value.

The wrapper for ?geqrf doesn't use this approach, it seems; it just defaults to 3n. Where does 3n come from? In the NAG manual, it looks like 64*n is used.

It is probably the absolute minimum value allowed.

[mdhaber]

Thanks @ilayn that's really helpful. I've noticed that for the other *_lwork routines most of the arguments are hidden, but some of the arguments are still (optionally) allowed. Can these affect the optimal lwork value? How do the authors of the wrappers know which to hide and which to keep as optional arguments?

Also, it looks like there is no ?geqrf_lwork even though ?geqrf is already wrapped. Should we create a ?geqrf_lwork routine while we're at it?

[ilayn]

That would be indeed nice. Often you can look at which arguments are passed to the lwork computation in the routine (ilaenv) and make them visible. For modern personal computers most of the times it shouldn't make any difference but better safe than sorry.

[mdhaber]

Finally, I'm curious - I played with scipy.linalg.dgeqrf with argument lwork=-1 and it gave me a value of lwork in work[0]. Then I called the routine with that value of lwork, and the routine exited (successfully) with a different value in work[0]. Supposedly, On exit, if INFO = 0, WORK(1) returns the optimal LWORK., so why did I get different values depending on whether I'm doing a workspace query or not?

[ilayn]

WORK array is the name LAPACK uses for "I use this array for scratch and temp values". In complex ones you also get iwork, rwork etc. And when you do LWORK=-1 it also uses the first entry of that array that to report the optimal value. That array size is optimal if you only give it big enough such that it uses Level 3 blocked algorithms comfortably and and small enough that you don't waste memory space.

That is the reason we wrap WORK as a scalar when we wrap <LAPACK_routine>_lwork so that it doesn't allocate 1000 array for just a scalar reporting.

TL;DR the value you see is often the intermediate value used for computing something except when LWORK=-1. You shouldn't see it because it should be hidden from user as it has no prescribed information unless explicitly stated by LAPACK devs.

[mdhaber]

I suppose I would have expected the documentation to say something like:
On exit, if INFO = 0 and LWORK = -1, WORK(1) returns the optimal LWORK.
But OK! Thanks for the tip about making work as a scalar in the *_lwork wrapper.

[Kai-Striega]

@mdhaber FYI I've pushed the geqrfp and geqrfp_lwork wrappers to my branch lapack_support_geqrfp.

A few (minor) points:

• ilaenv requires the problem dimensions to be specified, so the signature of the lwork function is geqrfp_lwork(m, n) not geqrfp_lwork(n) as you originally suggested.

• In many cases the value of lwork is calculated as n * nb, where n is the number of columns and nb is the blocksize. The blocksize is system dependent, on my system (Ubuntu/WSL x86_64) the blocksize is 32. This is much greater than the current default of 1 and I suspect the blocksize will be larger for most systems. Do you think it would be worth bumping to something larger?

[ilayn]

Minimum is defined in LAPACK docs. There are architectures that don't do blocked operations hence mininmum should be followed.

[mdhaber]

@swallan