<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def jacobi_prolongation_smoother(S, T, C, B, omega=4.0/3.0, degree=1, filter=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 : 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 -------- matrix([[ 1., 0.], [ 1., 0.], [ 1., 0.], [ 0., 1.], [ 0., 1.], [ 0., 1.]]) matrix([[ 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:
        # 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:
        # Carry out Jacobi, but after calculating the prolongator update, U,
        # apply satisfy constraints so that U*B = 0
        P = T
        for i 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 Sparsity 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 i in range(degree):
            P = P - (D_inv_S*P)

    return P