def sol0(M):
    return np.eye(M.shape[1]) * M[:, (np.newaxis), :]

def sol1(M):
    b = np.zeros((M.shape[0], M.shape[1], M.shape[1]))
diag = np.arange(M.shape[1])
b[:, (diag), (diag)] = M
return b
Numpy: Fastest way of computing diagonal for each row of a 2d array