def compute_BtBinv(B, C):
    """Create block inverses.

    Helper function that creates inv(B_i.T B_i) for each block row i in C,
    where B_i is B restricted to the sparsity pattern of block row i.

    Parameters
    ----------
    B : {array}
        (M,k) array, typically near-nullspace modes for coarse grid, i.e., B_c.
    C : {csr_matrix, bsr_matrix}
        Sparse NxM matrix, whose sparsity structure (i.e., matrix graph)
        is used to determine BtBinv.

    Returns
    -------
    BtBinv : {array}
        BtBinv[i] = inv(B_i.T B_i), where B_i is B restricted to the nonzero
        pattern of block row i in C.

    Examples
    --------
    >>> from numpy import array
    >>> from scipy.sparse import bsr_matrix
    >>> from pyamg.util.utils import compute_BtBinv
    >>> T = array([[ 1.,  0.],
    ...            [ 1.,  0.],
    ...            [ 0.,  .5],
    ...            [ 0.,  .25]])
    >>> T = bsr_matrix(T)
    >>> B = array([[1.],[2.]])
    >>> compute_BtBinv(B, T)
    array([[[ 1.  ]],
    <BLANKLINE>
           [[ 1.  ]],
    <BLANKLINE>
           [[ 0.25]],
    <BLANKLINE>
           [[ 0.25]]])

    Notes
    -----
    The principal calling routines are
    aggregation.smooth.energy_prolongation_smoother, and
    util.utils.filter_operator.

    BtBinv is used in the prolongation smoothing process that incorporates B
    into the span of prolongation with row-wise projection operators.  It is
    these projection operators that BtBinv is part of.

    """
    if not isspmatrix_bsr(C) and not isspmatrix_csr(C):
        raise TypeError('Expected bsr_matrix or csr_matrix for C')
    if C.shape[1] != B.shape[0]:
        raise TypeError('Expected matching dimensions such that C*B')

    # Problem parameters
    if isspmatrix_bsr(C):
        ColsPerBlock = C.blocksize[1]
        RowsPerBlock = C.blocksize[0]
    else:
        ColsPerBlock = 1
        RowsPerBlock = 1
    Ncoarse = C.shape[1]
    Nfine = C.shape[0]
    NullDim = B.shape[1]
    Nnodes = int(Nfine/RowsPerBlock)

    # Construct BtB
    BtBinv = np.zeros((Nnodes, NullDim, NullDim), dtype=B.dtype)
    BsqCols = sum(range(NullDim+1))
    Bsq = np.zeros((Ncoarse, BsqCols), dtype=B.dtype)
    counter = 0
    for i in range(NullDim):
        for j in range(i, NullDim):
            Bsq[:, counter] = np.conjugate(np.ravel(np.asarray(B[:, i]))) * \
                np.ravel(np.asarray(B[:, j]))
            counter = counter + 1
    # This specialized C-routine calculates (B.T B) for each row using Bsq
    pyamg.amg_core.calc_BtB(NullDim, Nnodes, ColsPerBlock,
                            np.ravel(np.asarray(Bsq)),
                            BsqCols, np.ravel(np.asarray(BtBinv)),
                            C.indptr, C.indices)

    # Invert each block of BtBinv, noting that amg_core.calc_BtB(...) returns
    # values in column-major form, thus necessitating the deep transpose
    #   This is the old call to a specialized routine, but lacks robustness
    #   pyamg.amg_core.pinv_array(np.ravel(BtBinv), Nnodes, NullDim, 'F')
    BtBinv = BtBinv.transpose((0, 2, 1)).copy()
    pinv_array(BtBinv)

    return BtBinv
Create block inverses.

    Helper function that creates inv(B_i.T B_i) for each block row i in C,
    where B_i is B restricted to the sparsity pattern of block row i.

    Parameters
    ----------
    B : {array}
        (M,k) array, typically near-nullspace modes for coarse grid, i.e., B_c.
    C : {csr_matrix, bsr_matrix}
        Sparse NxM matrix, whose sparsity structure (i.e., matrix graph)
        is used to determine BtBinv.

    Returns
    -------
    BtBinv : {array}
        BtBinv[i] = inv(B_i.T B_i), where B_i is B restricted to the nonzero
        pattern of block row i in C.

    Examples
    --------
    >>> from numpy import array
    >>> from scipy.sparse import bsr_matrix
    >>> from pyamg.util.utils import compute_BtBinv
    >>> T = array([[ 1.,  0.],
    ...            [ 1.,  0.],
    ...            [ 0.,  .5],
    ...            [ 0.,  .25]])
    >>> T = bsr_matrix(T)
    >>> B = array([[1.],[2.]])
    >>> compute_BtBinv(B, T)
    array([[[ 1.  ]],
    <BLANKLINE>
           [[ 1.  ]],
    <BLANKLINE>
           [[ 0.25]],
    <BLANKLINE>
           [[ 0.25]]])

    Notes
    -----
    The principal calling routines are
    aggregation.smooth.energy_prolongation_smoother, and
    util.utils.filter_operator.

    BtBinv is used in the prolongation smoothing process that incorporates B
    into the span of prolongation with row-wise projection operators.  It is
    these projection operators that BtBinv is part of.