def sym_bmat_mul(S, x, out):
    """
    Symmetric banded matrix times vector.

    Parameters
    ----------
    S : ndarray
        (b x n) array specifying symmetric banded matrix in
        "upper form".
    x : ndarray
        Vector of length n, multiplying S.
    out : ndarray
        Vector of length n, which is overwritten to store result.
    """

    b, n = S.shape

    for i in range(n):
        out[i] = S[-1, i] * x[i]

    for j in range(1, b):
        for i in range(n - j):
            out[i] += S[(-j - 1), j + i] * x[j + i]
            out[j + i] += S[(-j - 1), j + i] * x[i]

    return out