(B, C)
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.
Create block inverses.
Create block inverses .
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):
        cols_per_block = C.blocksize[1]
        rows_per_block = C.blocksize[0]
    else:
        cols_per_block = 1
        rows_per_block = 1
    Ncoarse = C.shape[1]
    Nfine = C.shape[0]
    NullDim = B.shape[1]
    Nnodes = int(Nfine/rows_per_block)

    # 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
    amg_core.calc_BtB(NullDim, Nnodes, cols_per_block,
                      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
    #   amg_core.pinv_array(np.ravel(BtBinv), Nnodes, NullDim, 'F')
    BtBinv = BtBinv.transpose((0, 2, 1)).copy()
    linalg.pinv_array(BtBinv)

    return BtBinv
def compute_BtBinv ( B , C ) : 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 ) : cols_per_block = C . blocksize [ 1 ] rows_per_block = C . blocksize [ 0 ] else : cols_per_block = 1 rows_per_block = 1 Ncoarse = C . shape [ 1 ] Nfine = C . shape [ 0 ] NullDim = B . shape [ 1 ] Nnodes = int ( Nfine / rows_per_block ) # 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 amg_core . calc_BtB ( NullDim , Nnodes , cols_per_block , 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 #   amg_core.pinv_array(np.ravel(BtBinv), Nnodes, NullDim, 'F') BtBinv = BtBinv . transpose ( ( 0 , 2 , 1 ) ) . copy ( ) linalg . pinv_array ( BtBinv ) return BtBinv