import torch


def resolvent(A, B, C, z):
    """
    A: (... N N)
    B: (... N)
    C: (... N)
    z: (... L)
    returns: (... L) represents C (z - A)^{-1} B
    """
    if A.is_complex() or B.is_complex() or C.is_complex() or z.is_complex():
        dtype = torch.cfloat
    else:
        dtype = torch.float
    N = A.shape[-1]
    I = torch.eye(N).to(device=A.device, dtype=dtype)
    A_ = I * z[...,None,None] - A.unsqueeze(-3)
    B = B.to(dtype)[..., None, :, None] # (... L N 1)
    r = torch.linalg.solve(A_, B).squeeze(-1) # (... L N)
    r = torch.sum(r*C.to(dtype).unsqueeze(-2), dim=-1)
    return r