(A, x, b, iterations=1, sweep='forward', omega=1.0,
                    Dinv=None)
Perform Gauss-Seidel iterations on the linear system A A.H x = b.

    Also known as Kaczmarz relaxation

    Parameters
    ----------
    A : csr_matrix
        Sparse NxN matrix
    x : ndarray
        Approximate solution (length N)
    b : ndarray
        Right-hand side (length N)
    iterations : int
        Number of iterations to perform
    sweep : {'forward','backward','symmetric'}
        Direction of sweep
    omega : float
        Relaxation parameter typically in (0, 2)
        if omega != 1.0, then algorithm becomes SOR on A A.H
    Dinv : ndarray
        Inverse of diag(A A.H),  (length N)

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

    References
    ----------
    .. [1] Brandt, Ta'asan.
       "Multigrid Method For Nearly Singular And Slightly Indefinite Problems."
       1985.  NASA Technical Report Numbers: ICASE-85-57; NAS 1.26:178026;
       NASA-CR-178026;

    .. [2] Kaczmarz.  Angenaeherte Aufloesung von Systemen Linearer
       Gleichungen. Bull. Acad.  Polon. Sci. Lett. A 35, 355-57.  1937

    Examples
    --------
    >>> # Use NE Gauss-Seidel as a Stand-Alone Solver
    >>> from pyamg.relaxation.relaxation import gauss_seidel_ne
    >>> from pyamg.gallery import poisson
    >>> from pyamg.util.linalg import norm
    >>> import numpy as np
    >>> A = poisson((10,10), format='csr')
    >>> x0 = np.zeros((A.shape[0],1))
    >>> b = np.ones((A.shape[0],1))
    >>> gauss_seidel_ne(A, x0, b, iterations=10, sweep='symmetric')
    >>> print(f'{norm(b-A*x0):2.4}')
    8.476
    >>> #
    >>> # Use NE Gauss-Seidel as the Multigrid Smoother
    >>> from pyamg import smoothed_aggregation_solver
    >>> sa = smoothed_aggregation_solver(A, B=np.ones((A.shape[0],1)),
    ...         coarse_solver='pinv', max_coarse=50,
    ...         presmoother=('gauss_seidel_ne', {'sweep' : 'symmetric'}),
    ...         postsmoother=('gauss_seidel_ne', {'sweep' : 'symmetric'}))
    >>> x0=np.zeros((A.shape[0],1))
    >>> residuals=[]
    >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Gauss-Seidel iterations on the linear system A A.H x = b.
Perform Gauss - Seidel iterations on the linear system A A . H x = b .
def gauss_seidel_ne(A, x, b, iterations=1, sweep='forward', omega=1.0,
                    Dinv=None):
    """Perform Gauss-Seidel iterations on the linear system A A.H x = b.

    Also known as Kaczmarz relaxation

    Parameters
    ----------
    A : csr_matrix
        Sparse NxN matrix
    x : ndarray
        Approximate solution (length N)
    b : ndarray
        Right-hand side (length N)
    iterations : int
        Number of iterations to perform
    sweep : {'forward','backward','symmetric'}
        Direction of sweep
    omega : float
        Relaxation parameter typically in (0, 2)
        if omega != 1.0, then algorithm becomes SOR on A A.H
    Dinv : ndarray
        Inverse of diag(A A.H),  (length N)

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

    References
    ----------
    .. [1] Brandt, Ta'asan.
       "Multigrid Method For Nearly Singular And Slightly Indefinite Problems."
       1985.  NASA Technical Report Numbers: ICASE-85-57; NAS 1.26:178026;
       NASA-CR-178026;

    .. [2] Kaczmarz.  Angenaeherte Aufloesung von Systemen Linearer
       Gleichungen. Bull. Acad.  Polon. Sci. Lett. A 35, 355-57.  1937

    Examples
    --------
    >>> # Use NE Gauss-Seidel as a Stand-Alone Solver
    >>> from pyamg.relaxation.relaxation import gauss_seidel_ne
    >>> from pyamg.gallery import poisson
    >>> from pyamg.util.linalg import norm
    >>> import numpy as np
    >>> A = poisson((10,10), format='csr')
    >>> x0 = np.zeros((A.shape[0],1))
    >>> b = np.ones((A.shape[0],1))
    >>> gauss_seidel_ne(A, x0, b, iterations=10, sweep='symmetric')
    >>> print(f'{norm(b-A*x0):2.4}')
    8.476
    >>> #
    >>> # Use NE Gauss-Seidel as the Multigrid Smoother
    >>> from pyamg import smoothed_aggregation_solver
    >>> sa = smoothed_aggregation_solver(A, B=np.ones((A.shape[0],1)),
    ...         coarse_solver='pinv', max_coarse=50,
    ...         presmoother=('gauss_seidel_ne', {'sweep' : 'symmetric'}),
    ...         postsmoother=('gauss_seidel_ne', {'sweep' : 'symmetric'}))
    >>> 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=['csr'])

    # Dinv for A*A.H
    if Dinv is None:
        Dinv = np.ravel(get_diagonal(A, norm_eq=2, inv=True))

    if sweep not in ('forward', 'backward', 'symmetric'):
        raise ValueError('valid sweep directions: "forward", "backward", and "symmetric"')

    if sweep == 'forward':
        row_start, row_stop, row_step = 0, len(x), 1
    elif sweep == 'backward':
        row_start, row_stop, row_step = len(x)-1, -1, -1
    elif sweep == 'symmetric':
        for _iter in range(iterations):
            gauss_seidel_ne(A, x, b, iterations=1, sweep='forward',
                            omega=omega, Dinv=Dinv)
            gauss_seidel_ne(A, x, b, iterations=1, sweep='backward',
                            omega=omega, Dinv=Dinv)
        return

    for _i in range(iterations):
        amg_core.gauss_seidel_ne(A.indptr, A.indices, A.data,
                                 x, b, row_start,
                                 row_stop, row_step, Dinv, omega)
def gauss_seidel_ne ( A , x , b , iterations = 1 , sweep = 'forward' , omega = 1.0 , Dinv = None ) : A , x , b = make_system ( A , x , b , formats = [ 'csr' ] ) # Dinv for A*A.H if Dinv is None : Dinv = np . ravel ( get_diagonal ( A , norm_eq = 2 , inv = True ) ) if sweep not in ( 'forward' , 'backward' , 'symmetric' ) : raise ValueError ( 'valid sweep directions: "forward", "backward", and "symmetric"' ) if sweep == 'forward' : row_start , row_stop , row_step = 0 , len ( x ) , 1 elif sweep == 'backward' : row_start , row_stop , row_step = len ( x ) - 1 , - 1 , - 1 elif sweep == 'symmetric' : for _iter in range ( iterations ) : gauss_seidel_ne ( A , x , b , iterations = 1 , sweep = 'forward' , omega = omega , Dinv = Dinv ) gauss_seidel_ne ( A , x , b , iterations = 1 , sweep = 'backward' , omega = omega , Dinv = Dinv ) return for _i in range ( iterations ) : amg_core . gauss_seidel_ne ( A . indptr , A . indices , A . data , x , b , row_start , row_stop , row_step , Dinv , omega )