def get_Cpt_params(A, Cnodes, AggOp, T):
    """Return C and F pts.

    Helper function that returns a dictionary of sparse matrices and arrays
    which allow us to easily operate on Cpts and Fpts separately.

    Parameters
    ----------
    A : {csr_matrix, bsr_matrix}
        Operator
    Cnodes : {array}
        Array of all root node indices.  This is an array of nodal indices,
        not degree-of-freedom indices.  If the blocksize of T is 1, then
        nodal indices and degree-of-freedom indices coincide.
    AggOp : {csr_matrix}
        Aggregation operator corresponding to A
    T : {bsr_matrix}
        Tentative prolongator based on AggOp

    Returns
    -------
    Dictionary containing these parameters:

    P_I : {bsr_matrix}
        Interpolation operator that carries out only simple injection from the
        coarse grid to fine grid Cpts nodes
    I_F : {bsr_matrix}
        Identity operator on Fpts, i.e., the action of this matrix zeros
        out entries in a vector at all Cpts, leaving Fpts untouched
    I_C : {bsr_matrix}
        Identity operator on Cpts nodes, i.e., the action of this matrix zeros
        out entries in a vector at all Fpts, leaving Cpts untouched
    Cpts : {array}
        An array of all root node dofs, corresponding to the F/C splitting
    Fpts : {array}
        An array of all non root node dofs, corresponding to the F/C splitting

    Examples
    --------
    >>> from numpy import array
    >>> from pyamg.util.utils import get_Cpt_params
    >>> from pyamg.gallery import poisson
    >>> from scipy.sparse import csr_matrix, bsr_matrix
    >>> A = poisson((10,), format='csr')
    >>> Cpts = array([3, 7])
    >>> AggOp = ([[ 1., 0.], [ 1., 0.],
    ...           [ 1., 0.], [ 1., 0.],
    ...           [ 1., 0.], [ 0., 1.],
    ...           [ 0., 1.], [ 0., 1.],
    ...           [ 0., 1.], [ 0., 1.]])
    >>> AggOp = csr_matrix(AggOp)
    >>> T = AggOp.copy().tobsr()
    >>> params = get_Cpt_params(A, Cpts, AggOp, T)
    >>> params['P_I'].todense()
    matrix([[ 0.,  0.],
            [ 0.,  0.],
            [ 0.,  0.],
            [ 1.,  0.],
            [ 0.,  0.],
            [ 0.,  0.],
            [ 0.,  0.],
            [ 0.,  1.],
            [ 0.,  0.],
            [ 0.,  0.]])

    Notes
    -----
    The principal calling routine is
    aggregation.smooth.energy_prolongation_smoother,
    which uses the Cpt_param dictionary for root-node style
    prolongation smoothing

    """
    if not isspmatrix_bsr(A) and not isspmatrix_csr(A):
        raise TypeError('Expected BSR or CSR matrix A')
    if not isspmatrix_csr(AggOp):
        raise TypeError('Expected CSR matrix AggOp')
    if not isspmatrix_bsr(T):
        raise TypeError('Expected BSR matrix T')
    if T.blocksize[0] != T.blocksize[1]:
        raise TypeError('Expected square blocksize for BSR matrix T')
    if A.shape[0] != A.shape[1]:
        raise TypeError('Expected square matrix A')
    if T.shape[0] != A.shape[0]:
        raise TypeError('Expected compatible dimensions for T and A,\
                         T.shape[0] = A.shape[0]')
    if Cnodes.shape[0] != AggOp.shape[1]:
        if AggOp.shape[1] > 1:
            raise TypeError('Number of columns in AggOp must equal number\
                             of Cnodes')

    if isspmatrix_bsr(A) and A.blocksize[0] > 1:
        # Expand the list of Cpt nodes to a list of Cpt dofs
        blocksize = A.blocksize[0]
        Cpts = np.repeat(blocksize*Cnodes, blocksize)
        for k in range(1, blocksize):
            Cpts[list(range(k, Cpts.shape[0], blocksize))] += k
    else:
        blocksize = 1
        Cpts = Cnodes
    Cpts = np.array(Cpts, dtype=int)

    # More input checking
    if Cpts.shape[0] != T.shape[1]:
        if T.shape[1] > blocksize:
            raise ValueError('Expected number of Cpts to match T.shape[1]')
    if blocksize != T.blocksize[0]:
        raise ValueError('Expected identical blocksize in A and T')
    if AggOp.shape[0] != int(T.shape[0]/blocksize):
        raise ValueError('Number of rows in AggOp must equal number of\
                          fine-grid nodes')

    # Create two maps, one for F points and one for C points
    ncoarse = T.shape[1]
    I_C = eye(A.shape[0], A.shape[1], format='csr')
    I_F = I_C.copy()
    I_F.data[Cpts] = 0.0
    I_F.eliminate_zeros()
    I_C = I_C - I_F
    I_C.eliminate_zeros()

    # Find Fpts, the complement of Cpts
    Fpts = I_F.indices.copy()

    # P_I only injects from Cpts on the coarse grid to the fine grid, but
    # because of it's later uses, it must have the CSC indices ordered as
    # in Cpts
    if I_C.nnz > 0:
        indices = Cpts.copy()
        indptr = np.arange(indices.shape[0]+1)
    else:
        indices = np.zeros((0,), dtype=T.indices.dtype)
        indptr = np.zeros((ncoarse+1,), dtype=T.indptr.dtype)

    P_I = csc_matrix((I_C.data.copy(), indices, indptr),
                     shape=(I_C.shape[0], ncoarse))
    P_I = P_I.tobsr(T.blocksize)

    # Use same blocksize as A
    if isspmatrix_bsr(A):
        I_C = I_C.tobsr(A.blocksize)
        I_F = I_F.tobsr(A.blocksize)
    else:
        I_C = I_C.tobsr(blocksize=(1, 1))
        I_F = I_F.tobsr(blocksize=(1, 1))

    return {'P_I': P_I, 'I_F': I_F, 'I_C': I_C, 'Cpts': Cpts, 'Fpts': Fpts}
Return C and F pts.

    Helper function that returns a dictionary of sparse matrices and arrays
    which allow us to easily operate on Cpts and Fpts separately.

    Parameters
    ----------
    A : {csr_matrix, bsr_matrix}
        Operator
    Cnodes : {array}
        Array of all root node indices.  This is an array of nodal indices,
        not degree-of-freedom indices.  If the blocksize of T is 1, then
        nodal indices and degree-of-freedom indices coincide.
    AggOp : {csr_matrix}
        Aggregation operator corresponding to A
    T : {bsr_matrix}
        Tentative prolongator based on AggOp

    Returns
    -------
    Dictionary containing these parameters:

    P_I : {bsr_matrix}
        Interpolation operator that carries out only simple injection from the
        coarse grid to fine grid Cpts nodes
    I_F : {bsr_matrix}
        Identity operator on Fpts, i.e., the action of this matrix zeros
        out entries in a vector at all Cpts, leaving Fpts untouched
    I_C : {bsr_matrix}
        Identity operator on Cpts nodes, i.e., the action of this matrix zeros
        out entries in a vector at all Fpts, leaving Cpts untouched
    Cpts : {array}
        An array of all root node dofs, corresponding to the F/C splitting
    Fpts : {array}
        An array of all non root node dofs, corresponding to the F/C splitting

    Examples
    --------
    >>> from numpy import array
    >>> from pyamg.util.utils import get_Cpt_params
    >>> from pyamg.gallery import poisson
    >>> from scipy.sparse import csr_matrix, bsr_matrix
    >>> A = poisson((10,), format='csr')
    >>> Cpts = array([3, 7])
    >>> AggOp = ([[ 1., 0.], [ 1., 0.],
    ...           [ 1., 0.], [ 1., 0.],
    ...           [ 1., 0.], [ 0., 1.],
    ...           [ 0., 1.], [ 0., 1.],
    ...           [ 0., 1.], [ 0., 1.]])
    >>> AggOp = csr_matrix(AggOp)
    >>> T = AggOp.copy().tobsr()
    >>> params = get_Cpt_params(A, Cpts, AggOp, T)
    >>> params['P_I'].todense()
    matrix([[ 0.,  0.],
            [ 0.,  0.],
            [ 0.,  0.],
            [ 1.,  0.],
            [ 0.,  0.],
            [ 0.,  0.],
            [ 0.,  0.],
            [ 0.,  1.],
            [ 0.,  0.],
            [ 0.,  0.]])

    Notes
    -----
    The principal calling routine is
    aggregation.smooth.energy_prolongation_smoother,
    which uses the Cpt_param dictionary for root-node style
    prolongation smoothing