def METHOD_NAME(A, product=None, modes=None, rtol=4e-8, atol=0., l2_err=0.):
    """SVD of a |VectorArray| using Gram-Schmidt orthogonalization.
    Viewing the |VectorArray| `A` as a `A.dim` x `len(A)` matrix, the
    return value of this method is the singular value decomposition of
    `A`, where the inner product on R^(`dim(A)`) is given by `product`
    and the inner product on R^(`len(A)`) is the Euclidean inner
    product.
    Parameters
    ----------
    A
        The |VectorArray| for which the SVD is to be computed.
    product
        Inner product |Operator| w.r.t. which the left singular vectors
        are computed.
    modes
        If not `None`, at most the first `modes` singular values and
        vectors are returned.
    rtol
        Singular values smaller than this value multiplied by the
        largest singular value are ignored.
    atol
        Singular values smaller than this value are ignored.
    l2_err
        Do not return more modes than needed to bound the
        l2-approximation error by this value. I.e. the number of
        returned modes is at most ::
            argmin_N { sum_{n=N+1}^{infty} s_n^2 <= l2_err^2 }
        where `s_n` denotes the n-th singular value.
    Returns
    -------
    U
        |VectorArray| of left singular vectors.
    s
        One-dimensional |NumPy array| of singular values.
    Vh
        |NumPy array| of right singular vectors.
    """
    assert isinstance(A, VectorArray)
    assert product is None or isinstance(product, Operator)
    if A.dim == 0 or len(A) == 0:
        return A.space.empty(), np.array([]), np.zeros((0, len(A)))
    logger = getLogger('pymor.algorithms.svd_va.qr_svd')
    with logger.block('Computing QR decomposition ...'):
        Q, R = gram_schmidt(A, product=product, return_R=True, check=False)
    with logger.block('Computing SVD of R ...'):
        U2, s, Vh = spla.svd(R, lapack_driver='gesvd')
    with logger.block('Choosing the number of modes ...'):
        tol = max(rtol * s[0], atol)
        above_tol = np.where(s >= tol)[0]
        if len(above_tol) == 0:
            return A.space.empty(), np.array([]), np.zeros((0, len(A)))
        last_above_tol = above_tol[-1]
        errs = np.concatenate((np.cumsum(s[::-1] ** 2)[::-1], [0.]))
        below_err = np.where(errs <= l2_err**2)[0]
        first_below_err = below_err[0]
        selected_modes = min(first_below_err, last_above_tol + 1)
        if modes is not None:
            selected_modes = min(selected_modes, modes)
        U2 = U2[:, :selected_modes]
        s = s[:selected_modes]
        Vh = Vh[:selected_modes]
    with logger.block(f'Computing left singular vectors ({selected_modes} modes) ...'):
        U = Q.lincomb(U2.T)
    return U, s, Vh