def cr(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, M=None,
       callback=None, residuals=None):
    """Conjugate Residual algorithm.

    Solves the linear system Ax = b. Left preconditioning is supported.
    The matrix A must be Hermitian symmetric (but not necessarily definite).

    Parameters
    ----------
    A : array, matrix, sparse matrix, LinearOperator
        n x n, linear system to solve
    b : array, matrix
        right hand side, shape is (n,) or (n,1)
    x0 : array, matrix
        initial guess, default is a vector of zeros
    tol : float
        relative convergence tolerance, i.e. tol is scaled by the
        preconditioner norm of r_0, or ||r_0||_M.
    maxiter : int
        maximum number of allowed iterations
    xtype : type
        dtype for the solution, default is automatic type detection
    M : array, matrix, sparse matrix, LinearOperator
        n x n, inverted preconditioner, i.e. solve M A x = M b.
    callback : function
        User-supplied function is called after each iteration as
        callback(xk), where xk is the current solution vector
    residuals : list
        residuals contains the residual norm history,
        including the initial residual.  The preconditioner norm
        is used, instead of the Euclidean norm.

    Returns
    -------
    (xNew, info)
    xNew : an updated guess to the solution of Ax = b
    info : halting status of cr

            ==  =======================================
            0   successful exit
            >0  convergence to tolerance not achieved,
                return iteration count instead.
            <0  numerical breakdown, or illegal input
            ==  =======================================

    Notes
    -----
    The LinearOperator class is in scipy.sparse.linalg.interface.
    Use this class if you prefer to define A or M as a mat-vec routine
    as opposed to explicitly constructing the matrix.  A.psolve(..) is
    still supported as a legacy.

    The 2-norm of the preconditioned residual is used both for halting and
    returned in the residuals list.

    Examples
    --------
    >>> from pyamg.krylov.cr import cr
    >>> from pyamg.util.linalg import norm
    >>> import numpy as np
    >>> from pyamg.gallery import poisson
    >>> A = poisson((10,10))
    >>> b = np.ones((A.shape[0],))
    >>> (x,flag) = cr(A,b, maxiter=2, tol=1e-8)
    >>> print norm(b - A*x)
    10.9370700187

    References
    ----------
    .. [1] Yousef Saad, "Iterative Methods for Sparse Linear Systems,
       Second Edition", SIAM, pp. 262-67, 2003
       http://www-users.cs.umn.edu/~saad/books.html

    """
    A, M, x, b, postprocess = make_system(A, M, x0, b)
    # n = len(b)
    # Ensure that warnings are always reissued from this function
    import warnings
    warnings.filterwarnings('always', module='pyamg\.krylov\._cr')

    # determine maxiter
    if maxiter is None:
        maxiter = int(1.3*len(b)) + 2
    elif maxiter < 1:
        raise ValueError('Number of iterations must be positive')

    # choose tolerance for numerically zero values
    # t = A.dtype.char
    # eps = np.finfo(np.float).eps
    # feps = np.finfo(np.single).eps
    # geps = np.finfo(np.longfloat).eps
    # _array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G': 2}
    # numerically_zero = {0: feps*1e3, 1: eps*1e6,
    #                     2: geps*1e6}[_array_precision[t]]

    # setup method
    r = b - A*x
    z = M*r
    p = z.copy()
    zz = np.inner(z.conjugate(), z)

    # use preconditioner norm
    normr = np.sqrt(zz)

    if residuals is not None:
        residuals[:] = [normr]  # initial residual

    # Check initial guess ( scaling by b, if b != 0,
    #   must account for case when norm(b) is very small)
    normb = norm(b)
    if normb == 0.0:
        normb = 1.0
    if normr < tol*normb:
        return (postprocess(x), 0)

    # Scale tol by ||r_0||_M
    if normr != 0.0:
        tol = tol*normr

    # How often should r be recomputed
    recompute_r = 8

    iter = 0

    Az = A*z
    rAz = np.inner(r.conjugate(), Az)
    Ap = A*p

    while True:

        rAz_old = rAz

        alpha = rAz / np.inner(Ap.conjugate(), Ap)       # 3
        x += alpha * p                           # 4

        if np.mod(iter, recompute_r) and iter > 0:       # 5
            r -= alpha * Ap
        else:
            r = b - A*x

        z = M*r

        Az = A*z
        rAz = np.inner(r.conjugate(), Az)

        beta = rAz/rAz_old                        # 6

        p *= beta                               # 7
        p += z

        Ap *= beta                               # 8
        Ap += Az

        iter += 1

        zz = np.inner(z.conjugate(), z)
        normr = np.sqrt(zz)                          # use preconditioner norm

        if residuals is not None:
            residuals.append(normr)

        if callback is not None:
            callback(x)

        if normr < tol:
            return (postprocess(x), 0)
        elif zz == 0.0:
            # important to test after testing normr < tol. rz == 0.0 is an
            # indicator of convergence when r = 0.0
            warn("\nSingular preconditioner detected in CR, ceasing \
                  iterations\n")
            return (postprocess(x), -1)

        if iter == maxiter:
            return (postprocess(x), iter)
def cr ( A , b , x0 = None , tol = 1e-5 , maxiter = None , xtype = None , M = None , callback = None , residuals = None ) : A , M , x , b , postprocess = make_system ( A , M , x0 , b ) import warnings warnings . filterwarnings ( 'always' , module = 'pyamg\.krylov\._cr' ) if maxiter is None : maxiter = int ( 1.3 * len ( b ) ) + 2 elif maxiter < 1 : raise ValueError ( 'Number of iterations must be positive' ) r = b - A * x z = M * r p = z . copy ( ) zz = np . inner ( z . conjugate ( ) , z ) normr = np . sqrt ( zz ) if residuals is not None : residuals [ : ] = [ normr ] normb = norm ( b ) if normb == 0.0 : normb = 1.0 if normr < tol * normb : return ( postprocess ( x ) , 0 ) if normr != 0.0 : tol = tol * normr recompute_r = 8 iter = 0 Az = A * z rAz = np . inner ( r . conjugate ( ) , Az ) Ap = A * p while True : rAz_old = rAz alpha = rAz / np . inner ( Ap . conjugate ( ) , Ap ) x += alpha * p if np . mod ( iter , recompute_r ) and iter > 0 : r -= alpha * Ap else : r = b - A * x z = M * r Az = A * z rAz = np . inner ( r . conjugate ( ) , Az ) beta = rAz / rAz_old p *= beta p += z Ap *= beta Ap += Az iter += 1 zz = np . inner ( z . conjugate ( ) , z ) normr = np . sqrt ( zz ) if residuals is not None : residuals . append ( normr ) if callback is not None : callback ( x ) if normr < tol : return ( postprocess ( x ) , 0 ) elif zz == 0.0 : warn ( "\nSingular preconditioner detected in CR, ceasing \
                  iterations\n" ) return ( postprocess ( x ) , - 1 ) if iter == maxiter : return ( postprocess ( x ) , iter )
Conjugate Residual algorithm.

    Solves the linear system Ax = b. Left preconditioning is supported.
    The matrix A must be Hermitian symmetric (but not necessarily definite).

    Parameters
    ----------
    A : array, matrix, sparse matrix, LinearOperator
        n x n, linear system to solve
    b : array, matrix
        right hand side, shape is (n,) or (n,1)
    x0 : array, matrix
        initial guess, default is a vector of zeros
    tol : float
        relative convergence tolerance, i.e. tol is scaled by the
        preconditioner norm of r_0, or ||r_0||_M.
    maxiter : int
        maximum number of allowed iterations
    xtype : type
        dtype for the solution, default is automatic type detection
    M : array, matrix, sparse matrix, LinearOperator
        n x n, inverted preconditioner, i.e. solve M A x = M b.
    callback : function
        User-supplied function is called after each iteration as
        callback(xk), where xk is the current solution vector
    residuals : list
        residuals contains the residual norm history,
        including the initial residual.  The preconditioner norm
        is used, instead of the Euclidean norm.

    Returns
    -------
    (xNew, info)
    xNew : an updated guess to the solution of Ax = b
    info : halting status of cr

            ==  =======================================
            0   successful exit
            >0  convergence to tolerance not achieved,
                return iteration count instead.
            <0  numerical breakdown, or illegal input
            ==  =======================================

    Notes
    -----
    The LinearOperator class is in scipy.sparse.linalg.interface.
    Use this class if you prefer to define A or M as a mat-vec routine
    as opposed to explicitly constructing the matrix.  A.psolve(..) is
    still supported as a legacy.

    The 2-norm of the preconditioned residual is used both for halting and
    returned in the residuals list.

    Examples
    --------
    >>> from pyamg.krylov.cr import cr
    >>> from pyamg.util.linalg import norm
    >>> import numpy as np
    >>> from pyamg.gallery import poisson
    >>> A = poisson((10,10))
    >>> b = np.ones((A.shape[0],))
    >>> (x,flag) = cr(A,b, maxiter=2, tol=1e-8)
    >>> print norm(b - A*x)
    10.9370700187

    References
    ----------
    .. [1] Yousef Saad, "Iterative Methods for Sparse Linear Systems,
       Second Edition", SIAM, pp. 262-67, 2003
       http://www-users.cs.umn.edu/~saad/books.html