def METHOD_NAME(mat_in, n_v, n_h):
    """
    Returns a 4D array of matrices, splitting mat_in into
    n_v * n_h square sub-arrays.
    """
    return [np.hsplit(x, n_h) for x in np.vsplit(mat_in, n_v)]