Source code for pyamg.krylov._cr

from numpy import inner, mod, sqrt
from scipy.sparse.linalg.isolve.utils import make_system
from pyamg.util.linalg import norm
from warnings import warn

__docformat__ = "restructuredtext en"

__all__ = ['cr']


[docs]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 >>> from pyamg.gallery import poisson >>> A = poisson((10,10)) >>> b = numpy.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, xtype=None) # 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 = numpy.finfo(numpy.float).eps # feps = numpy.finfo(numpy.single).eps # geps = numpy.finfo(numpy.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 = inner(z.conjugate(), z) # use preconditioner norm normr = 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 = inner(r.conjugate(), Az) Ap = A*p while True: rAz_old = rAz alpha = rAz / inner(Ap.conjugate(), Ap) # 3 x += alpha * p # 4 if mod(iter, recompute_r) and iter > 0: # 5 r -= alpha * Ap else: r = b - A*x z = M*r Az = A*z rAz = inner(r.conjugate(), Az) beta = rAz/rAz_old # 6 p *= beta # 7 p += z Ap *= beta # 8 Ap += Az iter += 1 zz = inner(z.conjugate(), z) normr = 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)
if __name__ == '__main__': # from numpy import diag # A = random((4,4)) # A = A*A.transpose() + diag([10,10,10,10]) # b = random((4,1)) # x0 = random((4,1)) from pyamg.gallery import stencil_grid from numpy.random import random import time from pyamg.krylov._gmres import gmres A = stencil_grid([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], (100, 100), dtype=float, format='csr') b = random((A.shape[0],)) x0 = random((A.shape[0],)) print '\n\nTesting CR with %d x %d 2D Laplace Matrix' % \ (A.shape[0], A.shape[0]) t1 = time.time() r = [] (x, flag) = cr(A, b, x0, tol=1e-8, maxiter=100, residuals=r) t2 = time.time() print '%s took %0.3f ms' % ('cr', (t2-t1)*1000.0) print 'norm = %g' % (norm(b - A*x)) print 'info flag = %d' % (flag) t1 = time.time() r2 = [] (x, flag) = gmres(A, b, x0, tol=1e-8, maxiter=100, residuals=r2) t2 = time.time() print '%s took %0.3f ms' % ('gmres', (t2-t1)*1000.0) print 'norm = %g' % (norm(b - A*x)) print 'info flag = %d' % (flag) # from scipy.sparse.linalg.isolve import cg as icg # t1=time.time() # (y,flag) = icg(A,b,x0,tol=1e-8,maxiter=100) # t2=time.time() # print '\n%s took %0.3f ms' % ('linalg cg', (t2-t1)*1000.0) # print 'norm = %g'%(norm(b - A*y)) # print 'info flag = %d'%(flag)