Gram-Schmidt orthogonalization with blas/lapack 08:24 PM

This is how to orthonormalize a set of vectors using BLAS and Lapack routines:

• I assume that you have the vectors stored in a matrix V.


• Calculate the overlap matrix S_ij = <v_i|v_j> with BLAS dsyrk (zherk for
  complex vectors).


• Calculate the Cholesky decomposition of S, this is to calculate a triangular matrix U such that S=U'*U. The decomposition is done by the Lapack function dpotrf.


• Your orthogonalized vectors will be W = inv(U)*V, but do not invert U and then multiply, as this is inefficient and possibly numerically unstable. Directly use the BLAS function dtrsm that multiplies a matrix in place by the inverse of a triangular matrix. Just what we need.

And that is all, just three subroutine calls.