Repository: scipy/scipy
Pull request: ENH: Tiled QR wrappers for scipy.linalg.lapack
As per https://github.com/scipy/scipy/issues/10323 I have begun working on the wrappers. I believe the wrappers and documentation are complete.

I still need to write some tests (but they do appear to function according to my brief, informal testing).

Closes #10323 

Problem statement:
Adding more low-level LAPACK wrappers
I would like to be able to make use of the following LAPACK routines from `scipy.linalg.lapack` (which currently have Cython wrappers provided in `scipy.linalg.cython_lapack`): `geqrt`, `gemqrt`, `tpqrt` and `tpmqrt` for tiled QR factorisation. I suggest the following signatures (with work and info as necessary):

```
?geqrt(nb, a, overwrite_a) -> (a, t)
?gemqrt(side, trans, v, t, c, overwrite_c) -> (c)
?tpqrt(l, nb, a, b, overwrite_a, overwrite_b) -> (a, b, t)
?tpmqrt(side, trans, v, t, a, b, overwrite_a, overwrite_b) -> (a, b)
```

Since QR on F-contiguous is equivalent to LQ on C-contiguous, in principle the LQ factorisation equivalents `gelqt`, `gemlqt`, `tplqt`, `tpmlqt` could be called on C-contiguous arrays to avoid a transposition at the high level. However, these routines do not have the Cython wrappers, which I assume is due to their relatively recent introduction to LAPACK.

I am happy to work on a PR for getting at least the four tiled QR routines going, but I will need some guidance with what's involved to do so.

Hints from issue discussion:
@AidanGG Thanks for taking the time. 

These routines came with Lapack 3.4 which we only recently assumed as the minimum version to support hence their absence.

For getting acquainted with the wrappers gou can check, say, #8965 on what is involved and which files are needed.

Please let me know if you need any more information. 
Hi @ilayn, that is very helpful. I will try my hand at writing a PR.

Commit message: Add missing f2py dependency

Patch to review:
diff --git a/THANKS.txt b/THANKS.txt
index e36e51797783..c499992ae1b1 100644
--- a/THANKS.txt
+++ b/THANKS.txt
@@ -209,6 +209,7 @@ William Conner DiPaolo for improvements to the Clarkson-Woodruff transform.
 Forrest Collman for adding multi-target dijsktra to scipy.sparse.csgraph
 Carlos Ramos Carreño for adding support for relational attributes in loadarff.
 Jason M. Manley for documentation fixes.
+Aidan Dang for block QR wrappers in scipy.linalg.lapack.
 
 Institutions
 ------------
diff --git a/scipy/linalg/flapack_other.pyf.src b/scipy/linalg/flapack_other.pyf.src
index 784426e18d07..69dc2093ea27 100644
--- a/scipy/linalg/flapack_other.pyf.src
+++ b/scipy/linalg/flapack_other.pyf.src
@@ -401,6 +401,116 @@ subroutine <prefix2c>unmqr(side,trans,m,n,k,a,lda,tau,c,ldc,work,lwork,info)
 
 end subroutine <prefix2c>unmqr
 
+subroutine <prefix>geqrt(m,n,nb,a,lda,t,ldt,work,info)
+    ! a,t,info = geqrt(nb,a,[overwrite_a=0])
+    !
+    ! Computes a QR factorization with block size nb of a general matrix A,
+    ! using the compact WY representation for Q. T stores the upper triangular
+    ! block reflectors.
+
+    callstatement (*f2py_func)(&m,&n,&nb,a,&lda,t,&ldt,work,&info)
+    callprotoargument int*,int*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*
+
+    integer intent(hide),depend(a):: m = shape(a,0)
+    integer intent(hide),depend(a):: n = shape(a,1)
+    integer intent(in),depend(m,n),check(MIN(m,n)>=nb&&nb>=1):: nb
+    <ftype> dimension(m,n),intent(in,out,copy):: a
+    integer intent(hide),depend(a):: lda = MAX(1,shape(a,0))
+    <ftype> dimension(nb,MIN(m,n)),intent(out):: t
+    integer intent(hide),depend(t):: ldt = MAX(1,shape(t,0))
+    <ftype> dimension(nb,n),intent(hide,cache):: work
+    integer intent(out):: info
+
+end subroutine <prefix>geqrt
+
+subroutine <prefix>gemqrt(side,trans,m,n,k,nb,v,ldv,t,ldt,c,ldc,work,info)
+    ! c,info = gemqrt(side,trans,v,t,c,[overwrite_c=0])
+    !
+    ! Multiplies a general matrix C by the orthogonal matrix Q defined by the
+    ! elementary reflectors stored in matrix V and the upper triangular block
+    ! reflectors in matrix T. C may be multiplied by Q, its transpose (for real
+    ! matrices), or its adjoint (for complex matrices) from the left or right.
+
+    callstatement (*f2py_func)(side,trans,&m,&n,&k,&nb,v,&ldv,t,&ldt,c,&ldc,work,&info)
+    callprotoargument char*,char*,int*,int*,int*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*
+
+    character intent(in),check(*side=='L'||*side=='R'):: side
+    character intent(in),check(*trans=='N'||*trans==<'T','T','C','C'>):: trans
+    integer intent(hide),depend(c):: m = shape(c,0)
+    integer intent(hide),depend(c):: n = shape(c,1)
+    integer intent(hide),depend(m,n,v),check((*side=='L'?m:n)>=k&&k>=0):: k = shape(v,1)
+    integer intent(hide),depend(k,t),check(k>=nb&&nb>=1):: nb = shape(t,0)
+    <ftype> dimension((*side=='L'?m:n),k),intent(in):: v
+    integer intent(hide),depend(v):: ldv = MAX(1,shape(v,0))
+    <ftype> dimension(nb,k),intent(in):: t
+    integer intent(hide),depend(t):: ldt = MAX(1,shape(t,0))
+    <ftype> dimension(m,n),intent(in,out,copy):: c
+    integer intent(hide),depend(c):: ldc = MAX(1,shape(c,0))
+    <ftype> dimension((*side=='L'?n:m)*nb),intent(hide,cache):: work
+    integer intent(out):: info
+
+end subroutine <prefix>gemqrt
+
+subroutine <prefix>tpqrt(m,n,l,nb,a,lda,b,ldb,t,ldt,work,info)
+    ! a,b,t,info = tpqrt(l,nb,a,b,[overwrite_a=0,overwrite_b=0])
+    !
+    ! Computes a QR factorization with block size nb of a triangular-pentagonal
+    ! matrix consisting of square upper triangular matrix A and pentagonal
+    ! matrix B, in the compact WY representation. L is the order of the
+    ! trapezoidal part of matrix B. T stores the upper triangular block
+    ! reflectors.
+
+    callstatement (*f2py_func)(&m,&n,&l,&nb,a,&lda,b,&ldb,t,&ldt,work,&info)
+    callprotoargument int*,int*,int*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*
+
+    integer intent(hide),depend(b):: m = shape(b,0)
+    integer intent(hide),depend(b):: n = shape(b,1)
+    integer intent(in),depend(m,n),check(MIN(m,n)>=l&&l>=0):: l
+    integer intent(in),depend(n),check(n>=nb&&nb>=1):: nb
+    <ftype> dimension(n,n),intent(in,out,copy):: a
+    integer intent(hide),depend(a):: lda = MAX(1,shape(a,0))
+    <ftype> dimension(m,n),intent(in,out,copy):: b
+    integer intent(hide),depend(b):: ldb = MAX(1,shape(b,0))
+    <ftype> dimension(nb,n),intent(out):: t
+    integer intent(hide),depend(t):: ldt = MAX(1,shape(t,0))
+    <ftype> dimension(nb,n),intent(hide,cache):: work
+    integer intent(out):: info
+
+end subroutine <prefix>tpqrt
+
+subroutine <prefix>tpmqrt(side,trans,m,n,k,l,nb,v,ldv,t,ldt,a,lda,b,ldb,work,info)
+    ! a,b,info = tpmqrt(side,trans,l,v,t,a,b,[overwrite_a=0,overwrite_b=0])
+    !
+    ! Multiplies a general matrix C by the orthogonal matrix Q defined by the
+    ! elementary reflectors stored in the pentagonal matrix V and the upper
+    ! triangular block reflectors in matrix T. L is the order of the trapezoidal
+    ! part of matrix V. Matrix C consists of blocks A and B, and may be
+    ! multiplied by Q, its transpose (for real matrices), or its adjoint (for
+    ! complex matrices) from the left or right.
+
+    callstatement (*f2py_func)(side,trans,&m,&n,&k,&l,&nb,v,&ldv,t,&ldt,a,&lda,b,&ldb,work,&info)
+    callprotoargument char*,char*,int*,int*,int*,int*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*
+
+    character intent(in),check(*side=='L'||*side=='R'):: side
+    character intent(in),check(*trans=='N'||*trans==<'T','T','C','C'>):: trans
+    integer intent(hide),depend(b):: m = shape(b,0)
+    integer intent(hide),depend(b):: n = shape(b,1)
+    integer intent(hide),depend(t):: k = shape(t,1)
+    integer intent(in),depend(k),check(k>=l&&l>=0):: l
+    integer intent(hide),depend(k,t),check(k>=nb&&nb>=1):: nb = shape(t,0)
+    <ftype> dimension((*side=='L'?m:n),k),intent(in):: v
+    integer intent(hide),depend(v):: ldv = MAX(1,shape(v,0))
+    <ftype> dimension(nb,k),intent(in):: t
+    integer intent(hide),depend(t):: ldt = MAX(1,shape(t,0))
+    <ftype> dimension((*side=='L'?k:m),(*side=='L'?n:k)),intent(in,out,copy):: a
+    integer intent(hide),depend(a):: lda = MAX(1,shape(a,0))
+    <ftype> dimension(m,n),intent(in,out,copy):: b
+    integer intent(hide),depend(b):: ldb = MAX(1,shape(b,0))
+    <ftype> dimension((*side=='L'?n:m)*nb),intent(hide,cache):: work
+    integer intent(out):: info
+
+end subroutine <prefix>tpmqrt
+
 subroutine <prefix><or,or,un,un>mrz(side,trans,m,n,k,l,a,lda,nt,tau,c,ldc,work,lwork,info)
     !
     ! ?OR/UNMRZ overwrites the general complex M-by-N matrix C with
@@ -419,7 +529,7 @@ subroutine <prefix><or,or,un,un>mrz(side,trans,m,n,k,l,a,lda,nt,tau,c,ldc,work,l
     !
     callstatement (*f2py_func)(side,trans,&m,&n,&k,&l,a,&lda,tau,c,&ldc,work,&lwork,&info)
     callprotoargument char*,char*,int*,int*,int*,int*,<ctype>*,int*,<ctype>*,<ctype>*,int*,<ctype>*,int*,int*
-    
+
     character optional,intent(in),check(*side=='L'||*side=='R'):: side = 'L'
     character optional,intent(in),check(*trans=='N'||*trans==<'T','T','C','C'>):: trans = 'N'
     integer intent(hide),depend(c),check(m>=0):: m=shape(c,0)
@@ -457,7 +567,7 @@ subroutine <prefix><or,or,un,un>mrz_lwork(side,trans,m,n,k,l,a,lda,tau,c,ldc,wor
     fortranname <prefix><or,or,un,un>mrz
     callstatement (*f2py_func)(side,trans,&m,&n,&k,&l,&a,&lda,&tau,&c,&ldc,&work,&lwork,&info)
     callprotoargument char*,char*,int*,int*,int*,int*,<ctype>*,int*,<ctype>*,<ctype>*,int*,<ctype>*,int*,int*
-    
+
     character optional,intent(in),check(*side=='L'||*side=='R'):: side = 'L'
     character optional,intent(in),check(*trans=='N'||*trans==<'T','T','C','C'>):: trans = 'N'
     integer intent(in):: m
@@ -472,7 +582,7 @@ subroutine <prefix><or,or,un,un>mrz_lwork(side,trans,m,n,k,l,a,lda,tau,c,ldc,wor
     <ftype> intent(out):: work
     integer intent(hide):: lwork=-1
     integer intent(out):: info
-    
+
 end subroutine <prefix><or,or,un,un>mrz_lwork
 
 subroutine <prefix2>orgrq(m,n,k,a,tau,work,lwork,info)
@@ -1060,7 +1170,7 @@ subroutine <prefix><s,s,h,h>frk(transr,uplo,trans,n,k,nt,ka,alpha,a,lda,beta,c)
     ! (S/D)SFRK - (C/Z)HFRK performs one of the Sym/Hermitian rank--k operations
     !
     !    C := alpha*A*A**H + beta*C,  or C := alpha*A**H*A + beta*C
-    ! 
+    !
     ! where alpha and beta are real scalars, C is an n--by--n Sym/Hermitian
     ! matrix and A is an n--by--k matrix in the first case and a k--by--n
     ! matrix in the second case.
@@ -1088,7 +1198,7 @@ subroutine <prefix>tpttf(transr,uplo,n,nt,ap,arf,info)
     ! copies a triangular matrix from the standard packed format (TP) to the
     ! rectangular full packed format (TF).
     !
-    
+
     callstatement (*f2py_func)(transr,uplo,&n,ap,arf,&info)
     callprotoargument char*,char*,int*,<ctype>*,<ctype>*,int*
 
@@ -1107,7 +1217,7 @@ subroutine <prefix>tpttr(uplo,n,nt,ap,a,lda,info)
     ! TPTTR copies a triangular matrix from the standard packed format (TP)
     ! to the standard full format (TR).
     !
-    
+
     callstatement (*f2py_func)(uplo,&n,ap,a,&lda,&info)
     callprotoargument char*,int*,<ctype>*,<ctype>*,int*,int*
 
@@ -1126,7 +1236,7 @@ subroutine <prefix>tfttp(transr,uplo,n,nt,ap,arf,info)
     ! copies a triangular matrix from the standard packed format (TP) to the
     ! rectangular full packed format (TF).
     !
-    
+
     callstatement (*f2py_func)(transr,uplo,&n,arf,ap,&info)
     callprotoargument char*,char*,int*,<ctype>*,<ctype>*,int*
 
@@ -1143,7 +1253,7 @@ end subroutine <prefix>tfttp
 subroutine <prefix>tfttr(transr,uplo,n,nt,arf,a,lda,info)
     !
     ! TFTTR copies a triangular matrix from the rectangular full packed
-    ! format (TF) to the standard full format (TR). 
+    ! format (TF) to the standard full format (TR).
     !
     callstatement (*f2py_func)(transr,uplo,&n,arf,a,&lda,&info)
     callprotoargument char*,char*,int*,<ctype>*,<ctype>*,int*,int*
@@ -1166,7 +1276,7 @@ subroutine <prefix>trttf(transr,uplo,n,a,lda,arf,info)
     !
     callstatement (*f2py_func)(transr,uplo,&n,a,&lda,arf,&info)
     callprotoargument char*,char*,int*,<ctype>*,int*,<ctype>*,int*
-    
+
     character optional,intent(in),check(*transr=='N'||*transr==<'T','T','C','C'>):: transr = 'N'
     character optional,intent(in),check(*uplo=='U'||*uplo=='L'):: uplo = 'U'
     integer intent(hide),depend(a):: n = shape(a,1)
@@ -1180,11 +1290,11 @@ end subroutine <prefix>trttf
 subroutine <prefix>trttp(uplo,n,a,lda,ap,info)
     !
     ! TRTTP copies a triangular matrix from the standard full format (TR) to
-    ! the standard packed format (TP). 
+    ! the standard packed format (TP).
     !
     callstatement (*f2py_func)(uplo,&n,a,&lda,ap,&info)
     callprotoargument char*,int*,<ctype>*,int*,<ctype>*,int*
-    
+
     character optional, intent(int),check(*uplo=='U'||*uplo=='L'):: uplo = 'U'
     integer intent(hide),depend(a):: n = shape(a,1)
     integer intent(hide),depend(a):: lda = MAX(shape(a,0),1)
@@ -1245,7 +1355,7 @@ subroutine <prefix>pftrf(transr,uplo,n,nt,a,info)
     character optional,intent(in),check(*uplo=='U'||*uplo=='L'):: uplo = 'U'
     integer intent(in), check(n>=0):: n
     <ftype> dimension(nt),intent(in,out,copy,out=achol):: a
-    integer intent(hide),depend(a,n),check(nt==(n*(n+1)/2)):: nt=shape(a,0)    
+    integer intent(hide),depend(a,n),check(nt==(n*(n+1)/2)):: nt=shape(a,0)
     integer intent(out):: info
 
 end subroutine <prefix>pftrf
@@ -1264,7 +1374,7 @@ subroutine <prefix>pftri(transr,uplo,n,nt,a,info)
     character optional,intent(in),check(*uplo=='U'||*uplo=='L'):: uplo = 'U'
     integer intent(in), check(n>=0):: n
     <ftype> dimension(nt),intent(in,out,copy,out=ainv):: a
-    integer intent(hide),depend(a,n),check(nt==(n*(n+1)/2)):: nt=shape(a,0)  
+    integer intent(hide),depend(a,n),check(nt==(n*(n+1)/2)):: nt=shape(a,0)
     integer intent(out):: info
 
 end subroutine <prefix>pftri
@@ -1303,7 +1413,7 @@ subroutine <prefix>tzrzf(m,n,a,lda,tau,work,lwork,info)
     !
     callstatement (*f2py_func)(&m,&n,a,&lda,tau,work,&lwork,&info)
     callprotoargument int*,int*,<ctype>*,int*,<ctype>*,<ctype>*,int*,int*
-    
+
     integer intent(hide),depend(a):: m=shape(a,0)
     integer intent(hide),depend(a):: n=shape(a,1)
     integer intent(hide),depend(a):: lda=MAX(shape(a,0),1)
@@ -1312,7 +1422,7 @@ subroutine <prefix>tzrzf(m,n,a,lda,tau,work,lwork,info)
     <ftype> dimension(MAX(lwork,1)),depend(lwork),intent(hide,cache):: work
     integer optional, intent(in),depend(m),check(lwork>=m):: lwork = MAX(m,1)
     integer intent(out):: info
-    
+
 end subroutine<prefix>tzrzf
 
 subroutine <prefix>tzrzf_lwork(m,n,a,lda,tau,work,lwork,info)
diff --git a/scipy/linalg/lapack.py b/scipy/linalg/lapack.py
index aed26797a175..21bf15c5afd3 100644
--- a/scipy/linalg/lapack.py
+++ b/scipy/linalg/lapack.py
@@ -597,6 +597,26 @@
    cunmqr
    zunmqr
 
+   sgeqrt
+   dgeqrt
+   cgeqrt
+   zgeqrt
+
+   sgemqrt
+   dgemqrt
+   cgemqrt
+   zgemqrt
+
+   stpqrt
+   dtpqrt
+   ctpqrt
+   ztpqrt
+
+   stpmqrt
+   dtpmqrt
+   ctpmqrt
+   ztpmqrt
+
    cunmrz
    zunmrz
 
diff --git a/scipy/linalg/tests/test_lapack.py b/scipy/linalg/tests/test_lapack.py
index 902aa12012c9..5e6d10533733 100644
--- a/scipy/linalg/tests/test_lapack.py
+++ b/scipy/linalg/tests/test_lapack.py
@@ -1276,3 +1276,133 @@ def test_syconv():
         ldu, ipiv, info = trf(A, lower=0)
         a, e, info = syconv(ldu, ipiv, lower=0)
         assert_allclose(triu(a, 1), triu(U[perm, :], 1), atol=tol, rtol=0.)
+
+class TestBlockedQR(object):
+    """
+    Tests for the blocked QR factorization, namely through geqrt, gemqrt, tpqrt
+    and tpmqr.
+    """
+
+    def test_geqrt_gemqrt(self):
+        seed(1234)
+        for ind, dtype in enumerate(DTYPES):
+            n = 20
+
+            if ind > 1:
+                A = (rand(n, n) + rand(n, n)*1j).astype(dtype)
+            else:
+                A = (rand(n, n)).astype(dtype)
+
+            tol = 100*np.spacing(dtype(1.0).real)
+            geqrt, gemqrt = get_lapack_funcs(('geqrt', 'gemqrt'), dtype=dtype)
+
+            a, t, info = geqrt(n, A)
+            assert(info == 0)
+
+            # Extract elementary reflectors from lower triangle, adding the main
+            # diagonal of ones.
+            v = np.tril(a, -1) + np.eye(n, dtype=dtype)
+            # Generate the block Householder transform I - VTV^H
+            Q = np.eye(n, dtype=dtype) - v @ t @ v.T.conj()
+            R = np.triu(a)
+
+            # Test columns of Q are orthogonal
+            assert_allclose(Q.T.conj() @ Q, np.eye(n, dtype=dtype), atol=tol,
+                            rtol=0.)
+            assert_allclose(Q @ R, A, atol=tol, rtol=0.)
+
+            if ind > 1:
+                C = (rand(n, n) + rand(n, n)*1j).astype(dtype)
+                transpose = 'C'
+            else:
+                C = (rand(n, n)).astype(dtype)
+                transpose = 'T'
+
+            for side in ('L', 'R'):
+                for trans in ('N', transpose):
+                    c, info = gemqrt(side, trans, a, t, C)
+                    assert(info == 0)
+
+                    if trans == transpose:
+                        q = Q.T.conj()
+                    else:
+                        q = Q
+
+                    if side == 'L':
+                        qC = q @ C
+                    else:
+                        qC = C @ q
+
+                    assert_allclose(c, qC, atol=tol, rtol=0.)
+
+    def test_tpqrt_tpmqrt(self):
+        seed(1234)
+        for ind, dtype in enumerate(DTYPES):
+            n = 20
+
+            if ind > 1:
+                A = (rand(n, n) + rand(n, n)*1j).astype(dtype)
+                B = (rand(n, n) + rand(n, n)*1j).astype(dtype)
+            else:
+                A = (rand(n, n)).astype(dtype)
+                B = (rand(n, n)).astype(dtype)
+
+            tol = 100*np.spacing(dtype(1.0).real)
+            tpqrt, tpmqrt = get_lapack_funcs(('tpqrt', 'tpmqrt'), dtype=dtype)
+
+            # Test for the range of pentagonal B, from square to upper
+            # triangular
+            for l in (0, n // 2, n):
+                a, b, t, info = tpqrt(l, n, A, B)
+                assert(info == 0)
+
+                # Check that lower triangular portion of A has not been modified
+                assert_equal(np.tril(a, -1), np.tril(A, -1))
+                # Check that elements not part of the pentagonal portion of B
+                # have not been modified.
+                assert_equal(np.tril(b, l - n - 1), np.tril(B, l - n - 1))
+
+                # Extract pentagonal portion of B
+                B_pent, b_pent = np.triu(B, l - n), np.triu(b, l - n)
+
+                # Generate elementary reflectors
+                v = np.concatenate((np.eye(n, dtype=dtype), b_pent))
+                # Generate the block Householder transform I - VTV^H
+                Q = np.eye(2 * n, dtype=dtype) - v @ t @ v.T.conj()
+                R = np.concatenate((np.triu(a), np.zeros_like(a)))
+
+                # Test columns of Q are orthogonal
+                assert_allclose(Q.T.conj() @ Q, np.eye(2 * n, dtype=dtype),
+                                atol=tol, rtol=0.)
+                assert_allclose(Q @ R, np.concatenate((np.triu(A), B_pent)),
+                                atol=tol, rtol=0.)
+
+                if ind > 1:
+                    C = (rand(n, n) + rand(n, n)*1j).astype(dtype)
+                    D = (rand(n, n) + rand(n, n)*1j).astype(dtype)
+                    transpose = 'C'
+                else:
+                    C = (rand(n, n)).astype(dtype)
+                    D = (rand(n, n)).astype(dtype)
+                    transpose = 'T'
+
+                for side in ('L', 'R'):
+                    for trans in ('N', transpose):
+                        c, d, info = tpmqrt(side, trans, l, b, t, C, D)
+                        assert(info == 0)
+
+                        if trans == transpose:
+                            q = Q.T.conj()
+                        else:
+                            q = Q
+
+                        if side == 'L':
+                            cd = np.concatenate((c, d), axis=0)
+                            CD = np.concatenate((C, D), axis=0)
+                            qCD = q @ CD
+                        else:
+                            cd = np.concatenate((c, d), axis=1)
+                            CD = np.concatenate((C, D), axis=1)
+                            qCD = CD @ q
+
+                        assert_allclose(cd, qCD, atol=tol, rtol=0.)


Code review comments:
diff_hunk: @@ -401,6 +401,116 @@ subroutine <prefix2c>unmqr(side,trans,m,n,k,a,lda,tau,c,ldc,work,lwork,info)
 
 end subroutine <prefix2c>unmqr
 
+subroutine <prefix>geqrt(m,n,nb,a,lda,t,ldt,work,info)
+    ! a,t,info = geqrt(nb,a,[overwrite_a=0])
+    !
+    ! Computes a QR factorization with block size nb of a general matrix A,
+    ! using the compact WY representation for Q. T stores the upper triangular
+    ! block reflectors.
+
+    callstatement (*f2py_func)(&m,&n,&nb,a,&lda,t,&ldt,work,&info)
+    callprotoargument int*,int*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*
+
+    integer intent(hide),depend(a):: m = shape(a,0)
+    integer intent(hide),depend(a):: n = shape(a,1)
+    integer intent(in),depend(m,n),check(MIN(m,n)>=nb&&nb>=1):: nb
+    <ftype> dimension(m,n),intent(in,out,copy):: a
+    integer intent(hide),depend(a):: lda = MAX(1,shape(a,0))
+    <ftype> dimension(nb,MIN(m,n)),intent(out):: t
+    integer intent(hide),depend(t):: ldt = MAX(1,shape(t,0))
+    <ftype> dimension(nb,n),intent(hide,cache):: work
+    integer intent(out):: info
+
+end subroutine <prefix>geqrt
+
+subroutine <prefix>gemqrt(side,trans,m,n,k,nb,v,ldv,t,ldt,c,ldc,work,info)
+    ! c,info = gemqrt(side,trans,v,t,c,[overwrite_c=0])
+    !
+    ! Multiplies a general matrix C by the orthogonal matrix Q defined by the
+    ! elementary reflectors stored in matrix V and the upper triangular block
+    ! reflectors in matrix T. C may be multiplied by Q, its transpose (for real
+    ! matrices), or its adjoint (for complex matrices) from the left or right.
+
+    callstatement (*f2py_func)(side,trans,&m,&n,&k,&nb,v,&ldv,t,&ldt,c,&ldc,work,&info)
+    callprotoargument char*,char*,int*,int*,int*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*,<ctype>*,int*
+
+    character intent(in),check(*side=='L'||*side=='R'):: side
path: scipy/linalg/flapack_other.pyf.src
text: @user1:
These chars best be optional with default values. I know we are not perfectly consistent but we are trying to fix them as we wrap new functions into these files. 

@user1:
Having said that we also need a simple test that gives wrong chars as arguments and we expect them to fail. 

@author:
Great, I'm assuming `trans` should also have the default value `'N'`?

@user1:
Yes indeed. And the side is L

@author:
Alright, done, with tests, in the latest commit.