(A, C, B, Bf, BtBinv=None)
return A
Filter the matrix A according to the matrix graph of C.

    Ensure that the new, filtered A satisfies:  A_new*B = Bf.

    A : {csr_matrix, bsr_matrix}
        n x m matrix to filter
    C : {csr_matrix, bsr_matrix}
        n x m matrix representing the couplings in A to keep
    B : {array}
        m x k array of near nullspace vectors
    Bf : {array}
        n x k array of near nullspace vectors to place in span(A)
    BtBinv : {None, array}
        3 dimensional array such that,
        BtBinv[i] = pinv(B_i.H Bi), and B_i is B restricted
        to the neighborhood (with respect to the matrix graph
        of C) of dof of i.  If None is passed in, this array is
        computed internally.

    Returns
    -------
    A : sparse matrix updated such that sparsity structure of A now matches
    that of C, and that the relationship A*B = Bf holds.

    Notes
    -----
    This procedure allows for certain important modes (i.e., Bf) to be placed
    in span(A) by way of row-wise l2-projections that enforce the relationship
    A*B = Bf.  This is useful for maintaining certain modes (e.g., the
    constant) in the span of prolongation.

    This routine is primarily used in
    pyamg.aggregation.smooth.energy_prolongation_smoother, where it is used to
    generate a suitable initial guess for the energy-minimization process, when
    root-node style SA is used.  Essentially, the tentative prolongator, T, is
    processed by this routine to produce fine-grid nullspace vectors when
    multiplying coarse-grid nullspace vectors, i.e., T*B = Bf.  This is
    possible for any arbitrary vectors B and Bf, so long as the sparsity
    structure of T is rich enough.

    When generating initial guesses for root-node style prolongation operators,
    this function is usually called before pyamg.uti.utils.scale_T

    Examples
    --------
    >>> from numpy import ones, array
    >>> from scipy.sparse import csr_matrix
    >>> from pyamg.util.utils import filter_operator
    >>> A = array([ [1.,1,1],[1,1,1],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
    >>> C = array([ [1.,1,0],[1,1,0],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
    >>> B = ones((3,1))
    >>> Bf = ones((6,1))
    >>> filter_operator(csr_matrix(A), csr_matrix(C), B, Bf).toarray()
    array([[0.5, 0.5, 0. ],
           [0.5, 0.5, 0. ],
           [0. , 1. , 0. ],
           [0. , 1. , 0. ],
           [0. , 0. , 1. ],
           [0. , 0. , 1. ]])
Filter the matrix A according to the matrix graph of C.
Filter the matrix A according to the matrix graph of C .
def filter_operator(A, C, B, Bf, BtBinv=None):
    """Filter the matrix A according to the matrix graph of C.

    Ensure that the new, filtered A satisfies:  A_new*B = Bf.

    A : {csr_matrix, bsr_matrix}
        n x m matrix to filter
    C : {csr_matrix, bsr_matrix}
        n x m matrix representing the couplings in A to keep
    B : {array}
        m x k array of near nullspace vectors
    Bf : {array}
        n x k array of near nullspace vectors to place in span(A)
    BtBinv : {None, array}
        3 dimensional array such that,
        BtBinv[i] = pinv(B_i.H Bi), and B_i is B restricted
        to the neighborhood (with respect to the matrix graph
        of C) of dof of i.  If None is passed in, this array is
        computed internally.

    Returns
    -------
    A : sparse matrix updated such that sparsity structure of A now matches
    that of C, and that the relationship A*B = Bf holds.

    Notes
    -----
    This procedure allows for certain important modes (i.e., Bf) to be placed
    in span(A) by way of row-wise l2-projections that enforce the relationship
    A*B = Bf.  This is useful for maintaining certain modes (e.g., the
    constant) in the span of prolongation.

    This routine is primarily used in
    pyamg.aggregation.smooth.energy_prolongation_smoother, where it is used to
    generate a suitable initial guess for the energy-minimization process, when
    root-node style SA is used.  Essentially, the tentative prolongator, T, is
    processed by this routine to produce fine-grid nullspace vectors when
    multiplying coarse-grid nullspace vectors, i.e., T*B = Bf.  This is
    possible for any arbitrary vectors B and Bf, so long as the sparsity
    structure of T is rich enough.

    When generating initial guesses for root-node style prolongation operators,
    this function is usually called before pyamg.uti.utils.scale_T

    Examples
    --------
    >>> from numpy import ones, array
    >>> from scipy.sparse import csr_matrix
    >>> from pyamg.util.utils import filter_operator
    >>> A = array([ [1.,1,1],[1,1,1],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
    >>> C = array([ [1.,1,0],[1,1,0],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
    >>> B = ones((3,1))
    >>> Bf = ones((6,1))
    >>> filter_operator(csr_matrix(A), csr_matrix(C), B, Bf).toarray()
    array([[0.5, 0.5, 0. ],
           [0.5, 0.5, 0. ],
           [0. , 1. , 0. ],
           [0. , 1. , 0. ],
           [0. , 0. , 1. ],
           [0. , 0. , 1. ]])

    """
    # First preprocess the parameters
    Nfine = A.shape[0]
    if A.shape[0] != C.shape[0]:
        raise ValueError('A and C must be the same size')
    if A.shape[1] != C.shape[1]:
        raise ValueError('A and C must be the same size')

    if isspmatrix_bsr(C):
        isBSR = True
        cols_per_block = C.blocksize[1]
        rows_per_block = C.blocksize[0]
        Nnodes = int(Nfine/rows_per_block)
        if not isspmatrix_bsr(A):
            raise ValueError('A and C must either both be CSR or BSR')

        if (cols_per_block != A.blocksize[1]) or (rows_per_block != A.blocksize[0]):
            raise ValueError('A and C must have same BSR blocksizes')

    elif isspmatrix_csr(C):
        isBSR = False
        cols_per_block = 1
        rows_per_block = 1
        Nnodes = int(Nfine/rows_per_block)
        if not isspmatrix_csr(A):
            raise ValueError('A and C must either both be CSR or BSR')
    else:
        raise ValueError('A and C must either both be CSR or BSR')

    if len(Bf.shape) == 1:
        Bf = Bf.reshape(-1, 1)
    if Bf.shape[0] != A.shape[0]:
        raise ValueError('A and Bf must have the same first dimension')

    if len(B.shape) == 1:
        B = B.reshape(-1, 1)
    if B.shape[0] != A.shape[1]:
        raise ValueError('A and B must have matching dimensions such\
                          that A*B is computable')

    if B.shape[1] != Bf.shape[1]:
        raise ValueError('B and Bf must have the same second\
                          dimension')

    NullDim = B.shape[1]

    if A.dtype == int:
        A.data = np.array(A.data, dtype=float)
    if B.dtype == int:
        B.data = np.array(B.data, dtype=float)
    if Bf.dtype == int:
        Bf.data = np.array(Bf.data, dtype=float)
    if (A.dtype != B.dtype) or (A.dtype != Bf.dtype):
        raise TypeError('A, B and Bf must of the same dtype')

    # First, preprocess some values for filtering.  Construct array of
    # inv(Bi'Bi), where Bi is B restricted to row i's sparsity pattern in
    # C. This array is used multiple times in satisfy_constraints(...).
    if BtBinv is None:
        BtBinv = compute_BtBinv(B, C)

    # Filter A according to C's matrix graph
    C = C.copy()
    C.data[:] = 1
    A = A.multiply(C)
    # add explicit zeros to A wherever C is nonzero, but A is zero
    A = A.tocoo()
    C = C.tocoo()
    A.data = np.hstack((np.zeros(C.data.shape, dtype=A.dtype), A.data))
    A.row = np.hstack((C.row, A.row))
    A.col = np.hstack((C.col, A.col))
    if isBSR:
        A = A.tobsr((rows_per_block, cols_per_block))
    else:
        A = A.tocsr()

    # Calculate difference between A*B and Bf
    diff = A*B - Bf

    # Right multiply each row i of A with
    # A_i <--- A_i - diff_i*inv(B_i.T B_i)*Bi.T
    # where A_i, and diff_i denote restriction to just row i, and B_i denotes
    # restriction to multiple rows corresponding to the the allowed nz's for
    # row i in A_i.  A_i also represents just the nonzeros for row i.
    amg_core.satisfy_constraints_helper(rows_per_block, cols_per_block,
                                        Nnodes, NullDim,
                                        np.conjugate(np.ravel(B)),
                                        np.ravel(diff),
                                        np.ravel(BtBinv), A.indptr,
                                        A.indices, np.ravel(A.data))

    A.eliminate_zeros()
    return A
def filter_operator ( A , C , B , Bf , BtBinv = None ) : # First preprocess the parameters Nfine = A . shape [ 0 ] if A . shape [ 0 ] != C . shape [ 0 ] : raise ValueError ( 'A and C must be the same size' ) if A . shape [ 1 ] != C . shape [ 1 ] : raise ValueError ( 'A and C must be the same size' ) if isspmatrix_bsr ( C ) : isBSR = True cols_per_block = C . blocksize [ 1 ] rows_per_block = C . blocksize [ 0 ] Nnodes = int ( Nfine / rows_per_block ) if not isspmatrix_bsr ( A ) : raise ValueError ( 'A and C must either both be CSR or BSR' ) if ( cols_per_block != A . blocksize [ 1 ] ) or ( rows_per_block != A . blocksize [ 0 ] ) : raise ValueError ( 'A and C must have same BSR blocksizes' ) elif isspmatrix_csr ( C ) : isBSR = False cols_per_block = 1 rows_per_block = 1 Nnodes = int ( Nfine / rows_per_block ) if not isspmatrix_csr ( A ) : raise ValueError ( 'A and C must either both be CSR or BSR' ) else : raise ValueError ( 'A and C must either both be CSR or BSR' ) if len ( Bf . shape ) == 1 : Bf = Bf . reshape ( - 1 , 1 ) if Bf . shape [ 0 ] != A . shape [ 0 ] : raise ValueError ( 'A and Bf must have the same first dimension' ) if len ( B . shape ) == 1 : B = B . reshape ( - 1 , 1 ) if B . shape [ 0 ] != A . shape [ 1 ] : raise ValueError ( 'A and B must have matching dimensions such\
                          that A*B is computable' ) if B . shape [ 1 ] != Bf . shape [ 1 ] : raise ValueError ( 'B and Bf must have the same second\
                          dimension' ) NullDim = B . shape [ 1 ] if A . dtype == int : A . data = np . array ( A . data , dtype = float ) if B . dtype == int : B . data = np . array ( B . data , dtype = float ) if Bf . dtype == int : Bf . data = np . array ( Bf . data , dtype = float ) if ( A . dtype != B . dtype ) or ( A . dtype != Bf . dtype ) : raise TypeError ( 'A, B and Bf must of the same dtype' ) # First, preprocess some values for filtering.  Construct array of # inv(Bi'Bi), where Bi is B restricted to row i's sparsity pattern in # C. This array is used multiple times in satisfy_constraints(...). if BtBinv is None : BtBinv = compute_BtBinv ( B , C ) # Filter A according to C's matrix graph C = C . copy ( ) C . data [ : ] = 1 A = A . multiply ( C ) # add explicit zeros to A wherever C is nonzero, but A is zero A = A . tocoo ( ) C = C . tocoo ( ) A . data = np . hstack ( ( np . zeros ( C . data . shape , dtype = A . dtype ) , A . data ) ) A . row = np . hstack ( ( C . row , A . row ) ) A . col = np . hstack ( ( C . col , A . col ) ) if isBSR : A = A . tobsr ( ( rows_per_block , cols_per_block ) ) else : A = A . tocsr ( ) # Calculate difference between A*B and Bf diff = A * B - Bf # Right multiply each row i of A with # A_i <--- A_i - diff_i*inv(B_i.T B_i)*Bi.T # where A_i, and diff_i denote restriction to just row i, and B_i denotes # restriction to multiple rows corresponding to the the allowed nz's for # row i in A_i.  A_i also represents just the nonzeros for row i. amg_core . satisfy_constraints_helper ( rows_per_block , cols_per_block , Nnodes , NullDim , np . conjugate ( np . ravel ( B ) ) , np . ravel ( diff ) , np . ravel ( BtBinv ) , A . indptr , A . indices , np . ravel ( A . data ) ) A . eliminate_zeros ( ) return A