linalg: QR factorization

[perazz]

Compute the QR factorization of a real or complex matrix: $A = Q R $, where q is orthonormal and r is upper-triangular. Matrix $A$ has size [m,n], with $m\ge n$. Based on LAPACK General QR factorization (*GEQRF) and ordered matrix output (*ORGQR, *UNGQR). Option for full or reduced factorization: given k = min(m,n), one can write $A = ( Q_1 Q_2 ) \cdot ( \frac{R_1}{0})$. The user may want the full problem or the reduced problem only $A = Q_1 R_1 $.

• [x] base implementation
• [x] tests
• [x] exclude unsupported xdp
• [x] documentation
• [x] submodule
• [x] examples
• [x] all pure subroutine interfaces
• [x] workspace size evaluation

Prior art

• Numpy: linalg.qr(a, mode={'reduced', 'complete', 'r', 'raw')
• Scipy: scipy.linalg.qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False, check_finite=True)

Proposed implementation

• call qr(A,Q,R [, overwrite_a] [, storage] [, err]) : pure subroutine interface
• call qr_space(A, lwork [, err]) query internal storage size for pre-allocation.

1) Special care was devoted to re-using internal storage such that a pure and totally allocation-less method is available, if the user provides pre-allocated working array storage. Because LAPACK internals require two steps ("raw" factorization + ordered matrix output), temporary storage is borrowed from either a, q, or r during the operations.

2) Both NumPy and SciPy have a cryptic UI that requires to provide a "method": full, reduced raw economic, r. I believe this is hard to understand. Instead, the current version proposes to decide it autonomously:

• If shape(Q)==[m,m] and shape(R)==[m,n], perform a full factorization, $A=QR$.
• If shape(Q)==[m,k] and shape(R)==[k,n], perform a "reduced" factorization, $A=Q_1R_1$.
• On insufficient/unfit size(s), raise an error.

In other words, the user decides what method is required based on the size of the arrays passed to qr.

cc: @fortran-lang/stdlib @jvdp1 @jalvesz @everythingfunctional

[perazz]

From Fortran Montly call:

• [x] require that Q and R have exact sizes for either the full or the reduced problem (do not just check that there is "enough" storage)