"""Iterative methods for solving linear systems""" from __future__ import division, print_function, absolute_import __all__ = ['bicg','bicgstab','cg','cgs','gmres','qmr'] from . import _iterative from scipy.sparse.linalg.interface import LinearOperator from scipy._lib.decorator import decorator from .utils import make_system from scipy._lib._util import _aligned_zeros from scipy._lib._threadsafety import non_reentrant _type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'} # Part of the docstring common to all iterative solvers common_doc1 = \ """ Parameters ---------- A : {sparse matrix, dense matrix, LinearOperator}""" common_doc2 = \ """b : {array, matrix} Right hand side of the linear system. Has shape (N,) or (N,1). Returns ------- x : {array, matrix} The converged solution. info : integer Provides convergence information: 0 : successful exit >0 : convergence to tolerance not achieved, number of iterations <0 : illegal input or breakdown Other Parameters ---------------- x0 : {array, matrix} Starting guess for the solution. tol : float Tolerance to achieve. The algorithm terminates when either the relative or the absolute residual is below `tol`. maxiter : integer Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved. M : {sparse matrix, dense matrix, LinearOperator} Preconditioner for A. The preconditioner should approximate the inverse of A. Effective preconditioning dramatically improves the rate of convergence, which implies that fewer iterations are needed to reach a given error tolerance. callback : function User-supplied function to call after each iteration. It is called as callback(xk), where xk is the current solution vector. xtype : {'f','d','F','D'} This parameter is deprecated -- avoid using it. The type of the result. If None, then it will be determined from A.dtype.char and b. If A does not have a typecode method then it will compute A.matvec(x0) to get a typecode. To save the extra computation when A does not have a typecode attribute use xtype=0 for the same type as b or use xtype='f','d','F',or 'D'. This parameter has been superseded by LinearOperator. """ def set_docstring(header, Ainfo, footer=''): def combine(fn): fn.__doc__ = '\n'.join((header, common_doc1, ' ' + Ainfo.replace('\n', '\n '), common_doc2, footer)) return fn return combine @set_docstring('Use BIConjugate Gradient iteration to solve A x = b', 'The real or complex N-by-N matrix of the linear system\n' 'It is required that the linear operator can produce\n' '``Ax`` and ``A^T x``.') @non_reentrant() def bicg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, M=None, callback=None): A,M,x,b,postprocess = make_system(A,M,x0,b,xtype) n = len(b) if maxiter is None: maxiter = n*10 matvec, rmatvec = A.matvec, A.rmatvec psolve, rpsolve = M.matvec, M.rmatvec ltr = _type_conv[x.dtype.char] revcom = getattr(_iterative, ltr + 'bicgrevcom') stoptest = getattr(_iterative, ltr + 'stoptest2') resid = tol ndx1 = 1 ndx2 = -1 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 work = _aligned_zeros(6*n,dtype=x.dtype) ijob = 1 info = 0 ftflag = True bnrm2 = -1.0 iter_ = maxiter while True: olditer = iter_ x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) if callback is not None and iter_ > olditer: callback(x) slice1 = slice(ndx1-1, ndx1-1+n) slice2 = slice(ndx2-1, ndx2-1+n) if (ijob == -1): if callback is not None: callback(x) break elif (ijob == 1): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(work[slice1]) elif (ijob == 2): work[slice2] *= sclr2 work[slice2] += sclr1*rmatvec(work[slice1]) elif (ijob == 3): work[slice1] = psolve(work[slice2]) elif (ijob == 4): work[slice1] = rpsolve(work[slice2]) elif (ijob == 5): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(x) elif (ijob == 6): if ftflag: info = -1 ftflag = False bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info) ijob = 2 if info > 0 and iter_ == maxiter and resid > tol: # info isn't set appropriately otherwise info = iter_ return postprocess(x), info @set_docstring('Use BIConjugate Gradient STABilized iteration to solve A x = b', 'The real or complex N-by-N matrix of the linear system\n' '``A`` must represent a hermitian, positive definite matrix') @non_reentrant() def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, M=None, callback=None): A,M,x,b,postprocess = make_system(A,M,x0,b,xtype) n = len(b) if maxiter is None: maxiter = n*10 matvec = A.matvec psolve = M.matvec ltr = _type_conv[x.dtype.char] revcom = getattr(_iterative, ltr + 'bicgstabrevcom') stoptest = getattr(_iterative, ltr + 'stoptest2') resid = tol ndx1 = 1 ndx2 = -1 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 work = _aligned_zeros(7*n,dtype=x.dtype) ijob = 1 info = 0 ftflag = True bnrm2 = -1.0 iter_ = maxiter while True: olditer = iter_ x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) if callback is not None and iter_ > olditer: callback(x) slice1 = slice(ndx1-1, ndx1-1+n) slice2 = slice(ndx2-1, ndx2-1+n) if (ijob == -1): if callback is not None: callback(x) break elif (ijob == 1): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(work[slice1]) elif (ijob == 2): work[slice1] = psolve(work[slice2]) elif (ijob == 3): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(x) elif (ijob == 4): if ftflag: info = -1 ftflag = False bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info) ijob = 2 if info > 0 and iter_ == maxiter and resid > tol: # info isn't set appropriately otherwise info = iter_ return postprocess(x), info @set_docstring('Use Conjugate Gradient iteration to solve A x = b', 'The real or complex N-by-N matrix of the linear system\n' '``A`` must represent a hermitian, positive definite matrix') @non_reentrant() def cg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, M=None, callback=None): A,M,x,b,postprocess = make_system(A,M,x0,b,xtype) n = len(b) if maxiter is None: maxiter = n*10 matvec = A.matvec psolve = M.matvec ltr = _type_conv[x.dtype.char] revcom = getattr(_iterative, ltr + 'cgrevcom') stoptest = getattr(_iterative, ltr + 'stoptest2') resid = tol ndx1 = 1 ndx2 = -1 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 work = _aligned_zeros(4*n,dtype=x.dtype) ijob = 1 info = 0 ftflag = True bnrm2 = -1.0 iter_ = maxiter while True: olditer = iter_ x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) if callback is not None and iter_ > olditer: callback(x) slice1 = slice(ndx1-1, ndx1-1+n) slice2 = slice(ndx2-1, ndx2-1+n) if (ijob == -1): if callback is not None: callback(x) break elif (ijob == 1): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(work[slice1]) elif (ijob == 2): work[slice1] = psolve(work[slice2]) elif (ijob == 3): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(x) elif (ijob == 4): if ftflag: info = -1 ftflag = False bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info) ijob = 2 if info > 0 and iter_ == maxiter and resid > tol: # info isn't set appropriately otherwise info = iter_ return postprocess(x), info @set_docstring('Use Conjugate Gradient Squared iteration to solve A x = b', 'The real-valued N-by-N matrix of the linear system') @non_reentrant() def cgs(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, M=None, callback=None): A,M,x,b,postprocess = make_system(A,M,x0,b,xtype) n = len(b) if maxiter is None: maxiter = n*10 matvec = A.matvec psolve = M.matvec ltr = _type_conv[x.dtype.char] revcom = getattr(_iterative, ltr + 'cgsrevcom') stoptest = getattr(_iterative, ltr + 'stoptest2') resid = tol ndx1 = 1 ndx2 = -1 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 work = _aligned_zeros(7*n,dtype=x.dtype) ijob = 1 info = 0 ftflag = True bnrm2 = -1.0 iter_ = maxiter while True: olditer = iter_ x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) if callback is not None and iter_ > olditer: callback(x) slice1 = slice(ndx1-1, ndx1-1+n) slice2 = slice(ndx2-1, ndx2-1+n) if (ijob == -1): if callback is not None: callback(x) break elif (ijob == 1): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(work[slice1]) elif (ijob == 2): work[slice1] = psolve(work[slice2]) elif (ijob == 3): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(x) elif (ijob == 4): if ftflag: info = -1 ftflag = False bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info) ijob = 2 if info > 0 and iter_ == maxiter and resid > tol: # info isn't set appropriately otherwise info = iter_ return postprocess(x), info @non_reentrant() def gmres(A, b, x0=None, tol=1e-5, restart=None, maxiter=None, xtype=None, M=None, callback=None, restrt=None): """ Use Generalized Minimal RESidual iteration to solve A x = b. Parameters ---------- A : {sparse matrix, dense matrix, LinearOperator} The real or complex N-by-N matrix of the linear system. b : {array, matrix} Right hand side of the linear system. Has shape (N,) or (N,1). Returns ------- x : {array, matrix} The converged solution. info : int Provides convergence information: * 0 : successful exit * >0 : convergence to tolerance not achieved, number of iterations * <0 : illegal input or breakdown Other parameters ---------------- x0 : {array, matrix} Starting guess for the solution (a vector of zeros by default). tol : float Tolerance to achieve. The algorithm terminates when either the relative or the absolute residual is below `tol`. restart : int, optional Number of iterations between restarts. Larger values increase iteration cost, but may be necessary for convergence. Default is 20. maxiter : int, optional Maximum number of iterations (restart cycles). Iteration will stop after maxiter steps even if the specified tolerance has not been achieved. xtype : {'f','d','F','D'} This parameter is DEPRECATED --- avoid using it. The type of the result. If None, then it will be determined from A.dtype.char and b. If A does not have a typecode method then it will compute A.matvec(x0) to get a typecode. To save the extra computation when A does not have a typecode attribute use xtype=0 for the same type as b or use xtype='f','d','F',or 'D'. This parameter has been superseded by LinearOperator. M : {sparse matrix, dense matrix, LinearOperator} Inverse of the preconditioner of A. M should approximate the inverse of A and be easy to solve for (see Notes). Effective preconditioning dramatically improves the rate of convergence, which implies that fewer iterations are needed to reach a given error tolerance. By default, no preconditioner is used. callback : function User-supplied function to call after each iteration. It is called as callback(rk), where rk is the current residual vector. restrt : int, optional DEPRECATED - use `restart` instead. See Also -------- LinearOperator Notes ----- A preconditioner, P, is chosen such that P is close to A but easy to solve for. The preconditioner parameter required by this routine is ``M = P^-1``. The inverse should preferably not be calculated explicitly. Rather, use the following template to produce M:: # Construct a linear operator that computes P^-1 * x. import scipy.sparse.linalg as spla M_x = lambda x: spla.spsolve(P, x) M = spla.LinearOperator((n, n), M_x) """ # Change 'restrt' keyword to 'restart' if restrt is None: restrt = restart elif restart is not None: raise ValueError("Cannot specify both restart and restrt keywords. " "Preferably use 'restart' only.") A,M,x,b,postprocess = make_system(A,M,x0,b,xtype) n = len(b) if maxiter is None: maxiter = n*10 if restrt is None: restrt = 20 restrt = min(restrt, n) matvec = A.matvec psolve = M.matvec ltr = _type_conv[x.dtype.char] revcom = getattr(_iterative, ltr + 'gmresrevcom') stoptest = getattr(_iterative, ltr + 'stoptest2') resid = tol ndx1 = 1 ndx2 = -1 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 work = _aligned_zeros((6+restrt)*n,dtype=x.dtype) work2 = _aligned_zeros((restrt+1)*(2*restrt+2),dtype=x.dtype) ijob = 1 info = 0 ftflag = True bnrm2 = -1.0 iter_ = maxiter old_ijob = ijob first_pass = True resid_ready = False iter_num = 1 while True: olditer = iter_ x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob) # if callback is not None and iter_ > olditer: # callback(x) slice1 = slice(ndx1-1, ndx1-1+n) slice2 = slice(ndx2-1, ndx2-1+n) if (ijob == -1): # gmres success, update last residual if resid_ready and callback is not None: callback(resid) resid_ready = False break elif (ijob == 1): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(x) elif (ijob == 2): work[slice1] = psolve(work[slice2]) if not first_pass and old_ijob == 3: resid_ready = True first_pass = False elif (ijob == 3): work[slice2] *= sclr2 work[slice2] += sclr1*matvec(work[slice1]) if resid_ready and callback is not None: callback(resid) resid_ready = False iter_num = iter_num+1 elif (ijob == 4): if ftflag: info = -1 ftflag = False bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info) old_ijob = ijob ijob = 2 if iter_num > maxiter: break if info >= 0 and resid > tol: # info isn't set appropriately otherwise info = maxiter return postprocess(x), info @non_reentrant() def qmr(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, M1=None, M2=None, callback=None): """Use Quasi-Minimal Residual iteration to solve A x = b Parameters ---------- A : {sparse matrix, dense matrix, LinearOperator} The real-valued N-by-N matrix of the linear system. It is required that the linear operator can produce ``Ax`` and ``A^T x``. b : {array, matrix} Right hand side of the linear system. Has shape (N,) or (N,1). Returns ------- x : {array, matrix} The converged solution. info : integer Provides convergence information: 0 : successful exit >0 : convergence to tolerance not achieved, number of iterations <0 : illegal input or breakdown Other Parameters ---------------- x0 : {array, matrix} Starting guess for the solution. tol : float Tolerance to achieve. The algorithm terminates when either the relative or the absolute residual is below `tol`. maxiter : integer Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved. M1 : {sparse matrix, dense matrix, LinearOperator} Left preconditioner for A. M2 : {sparse matrix, dense matrix, LinearOperator} Right preconditioner for A. Used together with the left preconditioner M1. The matrix M1*A*M2 should have better conditioned than A alone. callback : function User-supplied function to call after each iteration. It is called as callback(xk), where xk is the current solution vector. xtype : {'f','d','F','D'} This parameter is DEPRECATED -- avoid using it. The type of the result. If None, then it will be determined from A.dtype.char and b. If A does not have a typecode method then it will compute A.matvec(x0) to get a typecode. To save the extra computation when A does not have a typecode attribute use xtype=0 for the same type as b or use xtype='f','d','F',or 'D'. This parameter has been superseded by LinearOperator. See Also -------- LinearOperator """ A_ = A A,M,x,b,postprocess = make_system(A,None,x0,b,xtype) if M1 is None and M2 is None: if hasattr(A_,'psolve'): def left_psolve(b): return A_.psolve(b,'left') def right_psolve(b): return A_.psolve(b,'right') def left_rpsolve(b): return A_.rpsolve(b,'left') def right_rpsolve(b): return A_.rpsolve(b,'right') M1 = LinearOperator(A.shape, matvec=left_psolve, rmatvec=left_rpsolve) M2 = LinearOperator(A.shape, matvec=right_psolve, rmatvec=right_rpsolve) else: def id(b): return b M1 = LinearOperator(A.shape, matvec=id, rmatvec=id) M2 = LinearOperator(A.shape, matvec=id, rmatvec=id) n = len(b) if maxiter is None: maxiter = n*10 ltr = _type_conv[x.dtype.char] revcom = getattr(_iterative, ltr + 'qmrrevcom') stoptest = getattr(_iterative, ltr + 'stoptest2') resid = tol ndx1 = 1 ndx2 = -1 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 work = _aligned_zeros(11*n,x.dtype) ijob = 1 info = 0 ftflag = True bnrm2 = -1.0 iter_ = maxiter while True: olditer = iter_ x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) if callback is not None and iter_ > olditer: callback(x) slice1 = slice(ndx1-1, ndx1-1+n) slice2 = slice(ndx2-1, ndx2-1+n) if (ijob == -1): if callback is not None: callback(x) break elif (ijob == 1): work[slice2] *= sclr2 work[slice2] += sclr1*A.matvec(work[slice1]) elif (ijob == 2): work[slice2] *= sclr2 work[slice2] += sclr1*A.rmatvec(work[slice1]) elif (ijob == 3): work[slice1] = M1.matvec(work[slice2]) elif (ijob == 4): work[slice1] = M2.matvec(work[slice2]) elif (ijob == 5): work[slice1] = M1.rmatvec(work[slice2]) elif (ijob == 6): work[slice1] = M2.rmatvec(work[slice2]) elif (ijob == 7): work[slice2] *= sclr2 work[slice2] += sclr1*A.matvec(x) elif (ijob == 8): if ftflag: info = -1 ftflag = False bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info) ijob = 2 if info > 0 and iter_ == maxiter and resid > tol: # info isn't set appropriately otherwise info = iter_ return postprocess(x), info