I have a $p \times n$ matrix $B$ (where $n < p$) with orthonormal columns and would like to find a numerically efficient way to extend this matrix to get a complete $p$-dimensional orthonormal basis. In other words I want to compute a $p$ by $p-n$ complementary matrix $W$ so that
$[B|W]$
is orthonormal.
At the moment I am computing the QR factorisation of the rank $p-n$ projection matrix
$P=I-BB^T= QR$
using the LAPACK routines dgeqp3 and dorgqr and discarding the trailing $p$ columns of $Q$. I feel sure that there must be a more efficient method because in my application $B$ is almost a complete basis i.e $p-n$ is small.
I have wondered about building a complementary basis $R$ recursively by working through the columns of $P$. At each step adding a column $C$ to the current $R$ only if $\det([R|C]^T[R|C)\neq0$. After I have $p-n$ independent vectors in $R$ I could form $W$ as $W=(R^TR)^{-1/2} R$. However, I'm not sure if this is a numerically stable solution, or if there is another more efficient method.
Ideally I would like a method I can implement using standard LAPACK routines.