'Evaluate w = M^-1 v'

@staticmethod
def _solve(v, alpha, cs, ds):

    if (len(cs) == 0):
        return (v / alpha)
    (axpy, dotc) = get_blas_funcs(['axpy', 'dotc'], (cs[:1] + [v]))
    c0 = cs[0]
    A = (alpha * np.identity(len(cs), dtype=c0.dtype))
    for (i, d) in enumerate(ds):
        for (j, c) in enumerate(cs):
            A[(i, j)] += dotc(d, c)
    q = np.zeros(len(cs), dtype=c0.dtype)
    for (j, d) in enumerate(ds):
        q[j] = dotc(d, v)
    q /= alpha
    q = solve(A, q)
    w = (v / alpha)
    for (c, qc) in zip(cs, q):
        w = axpy(c, w, w.size, (- qc))
    return w
