While trying to find cell-centered gradients in finite volume method computation of incompressible fluid flow I get over-determined linear system. This is a well known "cell based least-square" gradient reconstruction.
I know in theory that by using least squares we search for a solution which minimizes error vector in 2-norm, and that we can solve the problem by calculating 'pseudo-inverse'
but I read that it is not typically how the problem is solved. Rather $QR$ decomposition is performed on original system matrix, and eventually we get
What I'm looking for is a fast algorithm for $QR$ decomposition of rectangular matrices that I can implement in Fortran, or any other suggestion related to particular problem that I mentioned.
System matrix is usually 6x3.