(A, x, b, coefficients, iterations=1)
Apply a polynomial smoother to the system Ax=b.

    Parameters
    ----------
    A : sparse matrix
        Sparse NxN matrix
    x : ndarray
        Approximate solution (length N)
    b : ndarray
        Right-hand side (length N)
    coefficients : array_like
        Coefficients of the polynomial.  See Notes section for details.
    iterations : int
        Number of iterations to perform

    Returns
    -------
    Nothing, x will be modified in place.

    Notes
    -----
    The smoother has the form  x[:] = x + p(A) (b - A*x) where p(A) is a
    polynomial in A whose scalar coefficients are specified (in descending
    order) by argument 'coefficients'.

    - Richardson iteration p(A) = c_0:
        polynomial_smoother(A, x, b, [c_0])

    - Linear smoother p(A) = c_1*A + c_0:
        polynomial_smoother(A, x, b, [c_1, c_0])

    - Quadratic smoother p(A) = c_2*A^2 + c_1*A + c_0:
        polynomial_smoother(A, x, b, [c_2, c_1, c_0])

    Here, Horner's Rule is applied to avoid computing A^k directly.

    For efficience, the method detects the case x = 0 one matrix-vector
    product is avoided (since (b - A*x) is b).

    Examples
    --------
    >>> # The polynomial smoother is not currently used directly
    >>> # in PyAMG.  It is only used by the chebyshev smoothing option,
    >>> # which automatically calculates the correct coefficients.
    >>> from pyamg.gallery import poisson
    >>> from pyamg.util.linalg import norm
    >>> import numpy as np
    >>> from pyamg.aggregation import smoothed_aggregation_solver
    >>> A = poisson((10,10), format='csr')
    >>> b = np.ones((A.shape[0],1))
    >>> sa = smoothed_aggregation_solver(A, B=np.ones((A.shape[0],1)),
    ...         coarse_solver='pinv', max_coarse=50,
    ...         presmoother=('chebyshev', {'degree':3, 'iterations':1}),
    ...         postsmoother=('chebyshev', {'degree':3, 'iterations':1}))
    >>> x0=np.zeros((A.shape[0],1))
    >>> residuals=[]
    >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Apply a polynomial smoother to the system Ax=b.
Apply a polynomial smoother to the system Ax = b .
def polynomial(A, x, b, coefficients, iterations=1):
    """Apply a polynomial smoother to the system Ax=b.

    Parameters
    ----------
    A : sparse matrix
        Sparse NxN matrix
    x : ndarray
        Approximate solution (length N)
    b : ndarray
        Right-hand side (length N)
    coefficients : array_like
        Coefficients of the polynomial.  See Notes section for details.
    iterations : int
        Number of iterations to perform

    Returns
    -------
    Nothing, x will be modified in place.

    Notes
    -----
    The smoother has the form  x[:] = x + p(A) (b - A*x) where p(A) is a
    polynomial in A whose scalar coefficients are specified (in descending
    order) by argument 'coefficients'.

    - Richardson iteration p(A) = c_0:
        polynomial_smoother(A, x, b, [c_0])

    - Linear smoother p(A) = c_1*A + c_0:
        polynomial_smoother(A, x, b, [c_1, c_0])

    - Quadratic smoother p(A) = c_2*A^2 + c_1*A + c_0:
        polynomial_smoother(A, x, b, [c_2, c_1, c_0])

    Here, Horner's Rule is applied to avoid computing A^k directly.

    For efficience, the method detects the case x = 0 one matrix-vector
    product is avoided (since (b - A*x) is b).

    Examples
    --------
    >>> # The polynomial smoother is not currently used directly
    >>> # in PyAMG.  It is only used by the chebyshev smoothing option,
    >>> # which automatically calculates the correct coefficients.
    >>> from pyamg.gallery import poisson
    >>> from pyamg.util.linalg import norm
    >>> import numpy as np
    >>> from pyamg.aggregation import smoothed_aggregation_solver
    >>> A = poisson((10,10), format='csr')
    >>> b = np.ones((A.shape[0],1))
    >>> sa = smoothed_aggregation_solver(A, B=np.ones((A.shape[0],1)),
    ...         coarse_solver='pinv', max_coarse=50,
    ...         presmoother=('chebyshev', {'degree':3, 'iterations':1}),
    ...         postsmoother=('chebyshev', {'degree':3, 'iterations':1}))
    >>> x0=np.zeros((A.shape[0],1))
    >>> residuals=[]
    >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)

    """
    A, x, b = make_system(A, x, b, formats=None)

    for _i in range(iterations):

        if norm(x) == 0:
            residual = b
        else:
            residual = (b - A*x)

        h = coefficients[0]*residual

        for c in coefficients[1:]:
            h = c*residual + A*h

        x += h
def polynomial ( A , x , b , coefficients , iterations = 1 ) : A , x , b = make_system ( A , x , b , formats = None ) for _i in range ( iterations ) : if norm ( x ) == 0 : residual = b else : residual = ( b - A * x ) h = coefficients [ 0 ] * residual for c in coefficients [ 1 : ] : h = c * residual + A * h x += h