Vectorize orthonormalization steps in gram_schmidt

[sdrave]

This refactors the gram_schmidt algorithm to orthonormalize all vectors to the right in a single vectorized operation. Depending on the use case, this can lead to significant speedup for VectorArrays like NumpyVectorArrays that support efficient vectorized operations.

[sdrave]

Some first quick check:

from IPython import get_ipython
ipython = get_ipython()


from pymor.basic import *
from pymor.algorithms.gram_schmidt import gram_schmidt_vec
U = NumpyVectorSpace(1000).random(100)
ipython.run_line_magic('timeit', 'gram_schmidt(U, reiterate=False)')
ipython.run_line_magic('timeit', 'gram_schmidt_vec(U, reiterate=False)')

reports on my Ryzen 6850U 133ms for the old implementation and 34ms for the new one. If I disable parallelization via OMP_NUM_THREADS=1 I get 76ms vs. 9ms. The results are the same up to machine precision.

[sdrave]

BTW, I think there is potential for quite a few optimizations on the Python side in case this turns out to be a limiting factor.

[maxbindhak]

While trying to add this to my tests i noticed that R gets cast into dtype 'O' (object) and after some debugging i could pinpoint it to those lines of code:

pymor/src/pymor/algorithms/gram_schmidt.py

Lines 218 to 221 in 87ac3e5

	p = A[i].inner(A[i+1:], product).ravel()
	A[i+1:].axpy(-p, A[i])
	common_dtype = np.promote_types(R.dtype, type(p))
	R = R.astype(common_dtype, copy=False)

Even thou the hack 87ac3e5 does not change the dtype of p, the common_dtype turns out to be 'O' and R gets cast into it.

For now i fixed it by simply casting R back into float64 at the end (i am only working with float64 anyway).

    if return_R:
        return A, R.astype(np.float64, copy=False)
    else:
        return A