(S, T, C, B, omega=4.0/3.0, degree=1,
                                 filter_entries=False, weighting='diagonal')
return P
Jacobi prolongation smoother.

    Parameters
    ----------
    S : csr_matrix, bsr_matrix
        Sparse NxN matrix used for smoothing.  Typically, A.
    T : csr_matrix, bsr_matrix
        Tentative prolongator
    C : csr_matrix, bsr_matrix
        Strength-of-connection matrix
    B : array
        Near nullspace modes for the coarse grid such that T*B
        exactly reproduces the fine grid near nullspace modes
    omega : scalar
        Damping parameter
    filter_entries : boolean
        If true, filter S before smoothing T.  This option can greatly control
        complexity.
    weighting : string
        'block', 'diagonal' or 'local' weighting for constructing the Jacobi D
        'local' Uses a local row-wise weight based on the Gershgorin estimate.
        Avoids any potential under-damping due to inaccurate spectral radius
        estimates.
        'block' uses a block diagonal inverse of A if A is BSR
        'diagonal' uses classic Jacobi with D = diagonal(A)

    Returns
    -------
    P : csr_matrix, bsr_matrix
        Smoothed (final) prolongator defined by P = (I - omega/rho(K) K) * T
        where K = diag(S)^-1 * S and rho(K) is an approximation to the
        spectral radius of K.

    Notes
    -----
    If weighting is not 'local', then results using Jacobi prolongation
    smoother are not precisely reproducible due to a random initial guess used
    for the spectral radius approximation.  For precise reproducibility,
    set numpy.random.seed(..) to the same value before each test.

    Examples
    --------
    >>> from pyamg.aggregation import jacobi_prolongation_smoother
    >>> from pyamg.gallery import poisson
    >>> from scipy.sparse import coo_matrix
    >>> import numpy as np
    >>> data = np.ones((6,))
    >>> row = np.arange(0,6)
    >>> col = np.kron([0,1],np.ones((3,)))
    >>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
    >>> T.toarray()
    array([[1., 0.],
           [1., 0.],
           [1., 0.],
           [0., 1.],
           [0., 1.],
           [0., 1.]])
    >>> A = poisson((6,),format='csr')
    >>> P = jacobi_prolongation_smoother(A,T,A,np.ones((2,1)))
    >>> P.toarray()
    array([[0.64930164, 0.        ],
           [1.        , 0.        ],
           [0.64930164, 0.35069836],
           [0.35069836, 0.64930164],
           [0.        , 1.        ],
           [0.        , 0.64930164]])
Jacobi prolongation smoother.
Jacobi prolongation smoother .
def jacobi_prolongation_smoother(S, T, C, B, omega=4.0/3.0, degree=1,
                                 filter_entries=False, weighting='diagonal'):
    """Jacobi prolongation smoother.

    Parameters
    ----------
    S : csr_matrix, bsr_matrix
        Sparse NxN matrix used for smoothing.  Typically, A.
    T : csr_matrix, bsr_matrix
        Tentative prolongator
    C : csr_matrix, bsr_matrix
        Strength-of-connection matrix
    B : array
        Near nullspace modes for the coarse grid such that T*B
        exactly reproduces the fine grid near nullspace modes
    omega : scalar
        Damping parameter
    filter_entries : boolean
        If true, filter S before smoothing T.  This option can greatly control
        complexity.
    weighting : string
        'block', 'diagonal' or 'local' weighting for constructing the Jacobi D
        'local' Uses a local row-wise weight based on the Gershgorin estimate.
        Avoids any potential under-damping due to inaccurate spectral radius
        estimates.
        'block' uses a block diagonal inverse of A if A is BSR
        'diagonal' uses classic Jacobi with D = diagonal(A)

    Returns
    -------
    P : csr_matrix, bsr_matrix
        Smoothed (final) prolongator defined by P = (I - omega/rho(K) K) * T
        where K = diag(S)^-1 * S and rho(K) is an approximation to the
        spectral radius of K.

    Notes
    -----
    If weighting is not 'local', then results using Jacobi prolongation
    smoother are not precisely reproducible due to a random initial guess used
    for the spectral radius approximation.  For precise reproducibility,
    set numpy.random.seed(..) to the same value before each test.

    Examples
    --------
    >>> from pyamg.aggregation import jacobi_prolongation_smoother
    >>> from pyamg.gallery import poisson
    >>> from scipy.sparse import coo_matrix
    >>> import numpy as np
    >>> data = np.ones((6,))
    >>> row = np.arange(0,6)
    >>> col = np.kron([0,1],np.ones((3,)))
    >>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
    >>> T.toarray()
    array([[1., 0.],
           [1., 0.],
           [1., 0.],
           [0., 1.],
           [0., 1.],
           [0., 1.]])
    >>> A = poisson((6,),format='csr')
    >>> P = jacobi_prolongation_smoother(A,T,A,np.ones((2,1)))
    >>> P.toarray()
    array([[0.64930164, 0.        ],
           [1.        , 0.        ],
           [0.64930164, 0.35069836],
           [0.35069836, 0.64930164],
           [0.        , 1.        ],
           [0.        , 0.64930164]])

    """
    # preprocess weighting
    if weighting == 'block':
        if sparse.isspmatrix_csr(S):
            weighting = 'diagonal'
        elif sparse.isspmatrix_bsr(S):
            if S.blocksize[0] == 1:
                weighting = 'diagonal'

    if filter_entries:
        # Implement filtered prolongation smoothing for the general case by
        # utilizing satisfy constraints

        if sparse.isspmatrix_bsr(S):
            numPDEs = S.blocksize[0]
        else:
            numPDEs = 1

        # Create a filtered S with entries dropped that aren't in C
        C = unamal(C, numPDEs, numPDEs)
        S = S.multiply(C)
        S.eliminate_zeros()

    if weighting == 'diagonal':
        # Use diagonal of S
        D_inv = get_diagonal(S, inv=True)
        D_inv_S = scale_rows(S, D_inv, copy=True)
        D_inv_S = (omega/approximate_spectral_radius(D_inv_S))*D_inv_S
    elif weighting == 'block':
        # Use block diagonal of S
        D_inv = get_block_diag(S, blocksize=S.blocksize[0], inv_flag=True)
        D_inv = sparse.bsr_matrix((D_inv, np.arange(D_inv.shape[0]),
                                   np.arange(D_inv.shape[0]+1)),
                                  shape=S.shape)
        D_inv_S = D_inv*S
        D_inv_S = (omega/approximate_spectral_radius(D_inv_S))*D_inv_S
    elif weighting == 'local':
        # Use the Gershgorin estimate as each row's weight, instead of a global
        # spectral radius estimate
        D = np.abs(S)*np.ones((S.shape[0], 1), dtype=S.dtype)
        D_inv = np.zeros_like(D)
        D_inv[D != 0] = 1.0 / np.abs(D[D != 0])

        D_inv_S = scale_rows(S, D_inv, copy=True)
        D_inv_S = omega*D_inv_S
    else:
        raise ValueError('Incorrect weighting option')

    if filter_entries:
        # Carry out Jacobi, but after calculating the prolongator update, U,
        # apply satisfy constraints so that U*B = 0
        P = T
        for _ in range(degree):
            U = (D_inv_S*P).tobsr(blocksize=P.blocksize)

            # Enforce U*B = 0 (1) Construct array of inv(Bi'Bi), where Bi is B
            # restricted to row i's sparsity pattern in pattern. This
            # array is used multiple times in satisfy_constraints(...).
            BtBinv = compute_BtBinv(B, U)

            # (2) Apply satisfy constraints
            satisfy_constraints(U, B, BtBinv)

            # Update P
            P = P - U

    else:
        # Carry out Jacobi as normal
        P = T
        for _ in range(degree):
            P = P - (D_inv_S*P)

    return P
def jacobi_prolongation_smoother ( S , T , C , B , omega = 4.0 / 3.0 , degree = 1 , filter_entries = False , weighting = 'diagonal' ) : # preprocess weighting if weighting == 'block' : if sparse . isspmatrix_csr ( S ) : weighting = 'diagonal' elif sparse . isspmatrix_bsr ( S ) : if S . blocksize [ 0 ] == 1 : weighting = 'diagonal' if filter_entries : # Implement filtered prolongation smoothing for the general case by # utilizing satisfy constraints if sparse . isspmatrix_bsr ( S ) : numPDEs = S . blocksize [ 0 ] else : numPDEs = 1 # Create a filtered S with entries dropped that aren't in C C = unamal ( C , numPDEs , numPDEs ) S = S . multiply ( C ) S . eliminate_zeros ( ) if weighting == 'diagonal' : # Use diagonal of S D_inv = get_diagonal ( S , inv = True ) D_inv_S = scale_rows ( S , D_inv , copy = True ) D_inv_S = ( omega / approximate_spectral_radius ( D_inv_S ) ) * D_inv_S elif weighting == 'block' : # Use block diagonal of S D_inv = get_block_diag ( S , blocksize = S . blocksize [ 0 ] , inv_flag = True ) D_inv = sparse . bsr_matrix ( ( D_inv , np . arange ( D_inv . shape [ 0 ] ) , np . arange ( D_inv . shape [ 0 ] + 1 ) ) , shape = S . shape ) D_inv_S = D_inv * S D_inv_S = ( omega / approximate_spectral_radius ( D_inv_S ) ) * D_inv_S elif weighting == 'local' : # Use the Gershgorin estimate as each row's weight, instead of a global # spectral radius estimate D = np . abs ( S ) * np . ones ( ( S . shape [ 0 ] , 1 ) , dtype = S . dtype ) D_inv = np . zeros_like ( D ) D_inv [ D != 0 ] = 1.0 / np . abs ( D [ D != 0 ] ) D_inv_S = scale_rows ( S , D_inv , copy = True ) D_inv_S = omega * D_inv_S else : raise ValueError ( 'Incorrect weighting option' ) if filter_entries : # Carry out Jacobi, but after calculating the prolongator update, U, # apply satisfy constraints so that U*B = 0 P = T for _ in range ( degree ) : U = ( D_inv_S * P ) . tobsr ( blocksize = P . blocksize ) # Enforce U*B = 0 (1) Construct array of inv(Bi'Bi), where Bi is B # restricted to row i's sparsity pattern in pattern. This # array is used multiple times in satisfy_constraints(...). BtBinv = compute_BtBinv ( B , U ) # (2) Apply satisfy constraints satisfy_constraints ( U , B , BtBinv ) # Update P P = P - U else : # Carry out Jacobi as normal P = T for _ in range ( degree ) : P = P - ( D_inv_S * P ) return P