Appending an existing QR decomposition: qrappend(q0, r0, Ap)

This algorithm computes the QR factorization of the matrix "A = q0*r0" with Ap appended at the end via the Modified Gram-Schmidt algorithm. It updates an existing QR decomposition without re-computation. The advantage over qrinsert is that it does not require that "Q must be square", that is, it allows the given QR decomposition to be of "economy size" / "thin form".
The input arguments q0 and r0 can be empty matrices, in that case the algorithm computes a new thin form QR of the matrix Ap.mM is the mass matrix defining a general inner product, setting it to "eye" or not supplying it recovers the regular dot product.

Example:

A = rand(10,5);
B = rand(10,2);
[q0, r0] = qr(A,0); % or qrappend([],[],A);
[q1, r1] = qrappend(q0,r0,B);