(A, b, x0=None, tol=1e-5, maxiter=400, return_solver=False,
          existing_solver=None, verb=True, residuals=None)
return x.reshape(b.shape)
Solve Ax=b.

    Solve the arbitrary system Ax=b with the best out-of-the box choice for a
    solver.  The matrix A can be non-Hermitian, indefinite, Hermitian
    positive-definite, complex, etc...  Generic and robust settings for
    smoothed_aggregation_solver(..) are used to invert A.


    Parameters
    ----------
    A : array, matrix, csr_matrix, bsr_matrix
        Matrix to invert, CSR or BSR format preferred for efficiency
    b : array
        Right hand side.
    x0 : array
        Initial guess (default random vector)
    tol : float
        Stopping criteria: relative residual r[k]/r[0] tolerance
    maxiter : int
        Stopping criteria: maximum number of allowable iterations
    return_solver : bool
        True: return the solver generated
    existing_solver : smoothed_aggregation_solver
        If instance of a multilevel solver, then existing_solver is used
        to invert A, thus saving time on setup cost.
    verb : bool
        If True, print verbose output during runtime
    residuals : list
        List to contain residual norms at each iteration.
        The preconditioned norm is used, namely
        ||r||_M = (M r, r)^(1/2) = (r, r)^(1/2)

    Returns
    -------
    x : array
        Solution to Ax = b
    ml : MultilevelSolver
        Optional return of the multilevel structure used for the solve

    Notes
    -----
    If calling solve(...) multiple times for the same matrix, A, solver reuse
    is easy and efficient.  Set "return_solver=True", and the return value will
    be a tuple, (x,ml), where ml is the solver used to invert A, and x is the
    solution to Ax=b.  Then, the next time solve(...) is called, set
    "existing_solver=ml".

    Examples
    --------
    >>> import numpy as np
    >>> from pyamg import solve
    >>> from pyamg.gallery import poisson
    >>> from pyamg.util.linalg import norm
    >>> A = poisson((40,40),format='csr')
    >>> b = np.array(np.arange(A.shape[0]), dtype=float)
    >>> x = solve(A,b,verb=False)
    >>> print(f'{norm(b - A*x)/norm(b):1.2e}')
    6.28e-06
Solve Ax=b.
Solve Ax = b .
def solve(A, b, x0=None, tol=1e-5, maxiter=400, return_solver=False,
          existing_solver=None, verb=True, residuals=None):
    """Solve Ax=b.

    Solve the arbitrary system Ax=b with the best out-of-the box choice for a
    solver.  The matrix A can be non-Hermitian, indefinite, Hermitian
    positive-definite, complex, etc...  Generic and robust settings for
    smoothed_aggregation_solver(..) are used to invert A.


    Parameters
    ----------
    A : array, matrix, csr_matrix, bsr_matrix
        Matrix to invert, CSR or BSR format preferred for efficiency
    b : array
        Right hand side.
    x0 : array
        Initial guess (default random vector)
    tol : float
        Stopping criteria: relative residual r[k]/r[0] tolerance
    maxiter : int
        Stopping criteria: maximum number of allowable iterations
    return_solver : bool
        True: return the solver generated
    existing_solver : smoothed_aggregation_solver
        If instance of a multilevel solver, then existing_solver is used
        to invert A, thus saving time on setup cost.
    verb : bool
        If True, print verbose output during runtime
    residuals : list
        List to contain residual norms at each iteration.
        The preconditioned norm is used, namely
        ||r||_M = (M r, r)^(1/2) = (r, r)^(1/2)

    Returns
    -------
    x : array
        Solution to Ax = b
    ml : MultilevelSolver
        Optional return of the multilevel structure used for the solve

    Notes
    -----
    If calling solve(...) multiple times for the same matrix, A, solver reuse
    is easy and efficient.  Set "return_solver=True", and the return value will
    be a tuple, (x,ml), where ml is the solver used to invert A, and x is the
    solution to Ax=b.  Then, the next time solve(...) is called, set
    "existing_solver=ml".

    Examples
    --------
    >>> import numpy as np
    >>> from pyamg import solve
    >>> from pyamg.gallery import poisson
    >>> from pyamg.util.linalg import norm
    >>> A = poisson((40,40),format='csr')
    >>> b = np.array(np.arange(A.shape[0]), dtype=float)
    >>> x = solve(A,b,verb=False)
    >>> print(f'{norm(b - A*x)/norm(b):1.2e}')
    6.28e-06

    """
    # Convert A to acceptable CSR/BSR format
    A = make_csr(A)

    # Generate solver if necessary
    if existing_solver is None:

        # Parameter dictionary for smoothed_aggregation_solver
        config = solver_configuration(A, B=None, verb=verb)
        # Generate solver
        existing_solver = solver(A, config)

    else:
        if existing_solver.levels[0].A.shape[0] != A.shape[0]:
            raise TypeError('Argument existing_solver must have level 0 matrix\
                             of same size as A')

    # Krylov acceleration depends on symmetry of A
    if existing_solver.levels[0].A.symmetry == 'hermitian':
        accel = 'cg'
    else:
        accel = 'gmres'

    # Initial guess
    if x0 is None:
        x0 = np.array(np.random.rand(A.shape[0],), dtype=A.dtype)

    # Callback function to print iteration number
    if verb:
        iteration = np.zeros((1,))
        print(f'    maxiter = {maxiter}')

        def callback(_x, iteration):
            iteration[0] = iteration[0] + 1
            print(f'    iteration {iteration[0]}')

        def callback2(x):
            return callback(x, iteration)
    else:
        callback2 = None

    # Solve with accelerated Krylov method
    x = existing_solver.solve(b, x0=x0, accel=accel, tol=tol, maxiter=maxiter,
                              callback=callback2, residuals=residuals)

    if verb:
        r0 = b - A * x0
        rk = b - A * x
        M = existing_solver.aspreconditioner()
        nr0 = np.sqrt(np.inner(np.conjugate(M * r0), r0))
        nrk = np.sqrt(np.inner(np.conjugate(M * rk), rk))
        print(f'  Residuals ||r_k||_M, ||r_0||_M = {nrk:1.2e}, {nr0:1.2e}')
        if np.abs(nr0) > 1e-15:
            ratio = nrk / nr0
            print(f'  Residual reduction ||r_k||_M/||r_0||_M = {ratio:1.2e}')

    if return_solver:
        return (x.reshape(b.shape), existing_solver)
    return x.reshape(b.shape)
def solve ( A , b , x0 = None , tol = 1e-5 , maxiter = 400 , return_solver = False , existing_solver = None , verb = True , residuals = None ) : # Convert A to acceptable CSR/BSR format A = make_csr ( A ) # Generate solver if necessary if existing_solver is None : # Parameter dictionary for smoothed_aggregation_solver config = solver_configuration ( A , B = None , verb = verb ) # Generate solver existing_solver = solver ( A , config ) else : if existing_solver . levels [ 0 ] . A . shape [ 0 ] != A . shape [ 0 ] : raise TypeError ( 'Argument existing_solver must have level 0 matrix\
                             of same size as A' ) # Krylov acceleration depends on symmetry of A if existing_solver . levels [ 0 ] . A . symmetry == 'hermitian' : accel = 'cg' else : accel = 'gmres' # Initial guess if x0 is None : x0 = np . array ( np . random . rand ( A . shape [ 0 ] , ) , dtype = A . dtype ) # Callback function to print iteration number if verb : iteration = np . zeros ( ( 1 , ) ) print ( f'    maxiter = {maxiter}' ) def callback ( _x , iteration ) : iteration [ 0 ] = iteration [ 0 ] + 1 print ( f'    iteration {iteration[0]}' ) def callback2 ( x ) : return callback ( x , iteration ) else : callback2 = None # Solve with accelerated Krylov method x = existing_solver . solve ( b , x0 = x0 , accel = accel , tol = tol , maxiter = maxiter , callback = callback2 , residuals = residuals ) if verb : r0 = b - A * x0 rk = b - A * x M = existing_solver . aspreconditioner ( ) nr0 = np . sqrt ( np . inner ( np . conjugate ( M * r0 ) , r0 ) ) nrk = np . sqrt ( np . inner ( np . conjugate ( M * rk ) , rk ) ) print ( f'  Residuals ||r_k||_M, ||r_0||_M = {nrk:1.2e}, {nr0:1.2e}' ) if np . abs ( nr0 ) > 1e-15 : ratio = nrk / nr0 print ( f'  Residual reduction ||r_k||_M/||r_0||_M = {ratio:1.2e}' ) if return_solver : return ( x . reshape ( b . shape ) , existing_solver ) return x . reshape ( b . shape )