@staticmethod<EOL><INDENT>def _solve(v, alpha, cs, ds):<DEDENT>
if len(cs) == <NUM_LIT:0>:<EOL><INDENT>return v/alpha<EOL><DEDENT>axpy, dotc = get_blas_funcs(['<STR_LIT>', '<STR_LIT>'], cs[:<NUM_LIT:1>] + [v])<EOL>c0 = cs[<NUM_LIT:0>]<EOL>A = alpha * np.identity(len(cs), dtype=c0.dtype)<EOL>for i, d in enumerate(ds):<EOL><INDENT>for j, c in enumerate(cs):<EOL><INDENT>A[i,j] += dotc(d, c)<EOL><DEDENT><DEDENT>q = np.zeros(len(cs), dtype=c0.dtype)<EOL>for j, d in enumerate(ds):<EOL><INDENT>q[j] = dotc(d, v)<EOL><DEDENT>q /= alpha<EOL>q = solve(A, q)<EOL>w = v/alpha<EOL>for c, qc in zip(cs, q):<EOL><INDENT>w = axpy(c, w, w.size, -qc)<EOL><DEDENT>return w<EOL>
Evaluate w = M^-1 v