3 Odl @s dZddddddgZddlZddlZd d lmZdd lmZd d l m Z dd l m Z ddl mZdddddZdZdZddZddZd/ddZeddd d!ed0d#dZed$d%ed1d&dZed'd(ed2d)dZed*d+ed3d,dZed4d-dZed5d.dZdS)6z,Iterative methods for solving linear systemsbicgbicgstabcgcgsgmresqmrN) _iterative)LinearOperator) make_system)_aligned_zeros) non_reentrantsdcz)frFDzH Parameters ---------- A : {sparse matrix, dense matrix, LinearOperator}awb : {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, atol : float, optional Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``. The default for ``atol`` is ``'legacy'``, which emulates a different legacy behavior. .. warning:: The default value for `atol` will be changed in a future release. For future compatibility, specify `atol` explicitly. 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. cCs(tjj|}||kr|dfS|dfSdS)z; Successful termination condition for the solvers. rrN)nplinalgnorm)ZresidualatolresidrX/var/www/html/virt/lib64/python3.6/site-packages/scipy/sparse/linalg/isolve/iterative.py _stoptestCs rcCsz|dkr$tjdj|dtddd}t|}|dkr`|}||krFdS|dkrR|S|t|Sntt||t|SdS) a Parse arguments for absolute tolerance in termination condition. Parameters ---------- tol, atol : object The arguments passed into the solver routine by user. bnrm2 : float 2-norm of the rhs vector. get_residual : callable Callable ``get_residual()`` that returns the initial value of the residual. routine_name : str Name of the routine. Na scipy.sparse.linalg.{name} called without specifying `atol`. The default value will be changed in a future release. For compatibility, specify a value for `atol` explicitly, e.g., ``{name}(..., atol=0)``, or to retain the old behavior ``{name}(..., atol='legacy')``)name)category stacklevellegacyexitr)warningswarnformatDeprecationWarningfloatmax)tolrbnrm2 get_residualZ routine_namerrrr _get_atolNs r,0csfdd}|S)Ncs&djtdjddtf|_|S)N z z )join common_doc1replace common_doc2__doc__)fn)Ainfofooterheaderrrcombinexs zset_docstring..combiner)r8r6r7Z atol_defaultr9r)r6r7r8r set_docstringwsr:z7Use BIConjugate Gradient iteration to solve ``Ax = b``.zThe real or complex N-by-N matrix of the linear system. Alternatively, ``A`` can be a linear operator which can produce ``Ax`` and ``A^T x`` using, e.g., ``scipy.sparse.linalg.LinearOperator``.a+ Examples -------- >>> from scipy.sparse import csc_matrix >>> from scipy.sparse.linalg import bicg >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float) >>> b = np.array([2, 4, -1], dtype=float) >>> x, exitCode = bicg(A, b) >>> print(exitCode) # 0 indicates successful convergence 0 >>> np.allclose(A.dot(x), b) True )r7h㈵>c st|||\}}}t} |dkr0| d}|j|j} |j|j} } tjj} tt| d}fdd}t ||t j j |d}|dkr|dfS|}d}d}t d | jd }d}d}d }|}x|}|||||||| \ }}}}}}}}|dk r||kr|t|d|d| }t|d|d| }|dkrl|dk rf|Pn|dkr|||9<|||||7<n|d kr|||9<|||| ||7<n|d kr| ||||<nz|dkr| ||||<n^|dkrH|||9<|||7<n*|d krr|r`d}d}t|||\}}d }qW|dkr||kr||k r|}||fS)N Z bicgrevcomcstjjS)N)rrrr)bmatvecxrrszbicg..rr"rr)dtypeTrFrFrF)r lenr>rmatvec _type_convrBchargetattrr r,rrrr slicer)Ar=x0r)maxiterMcallbackr postprocessnrHpsolverpsolveltrrevcomr+rndx1ndx2workijobinfoftflagiter_olditersclr1sclr2slice1slice2r)r=r>r?rrsj  *         zBUse BIConjugate Gradient STABilized iteration to solve ``Ax = b``.zThe real or complex N-by-N matrix of the linear system. Alternatively, ``A`` can be a linear operator which can produce ``Ax`` using, e.g., ``scipy.sparse.linalg.LinearOperator``.c sDt|||\}}}t} |dkr0| d}|j|j} tjj} tt| d} fdd} t||t j j | d}|dkr|dfS|}d}d}t d | jd }d}d}d }|}xP|}| ||||||| \ }}}}}}}}|dk r||kr|t |d|d| }t |d|d| }|dkrZ|dk rV|Pn|dkr|||9<|||||7<nz|d kr| ||||<n^|d kr|||9<|||7<n*|dkr |rd}d}t|||\}}d }qW|dkr8||kr8||k r8|}||fS)Nr<ZbicgstabrevcomcstjjS)N)rrrr)r=r>r?rrr@szbicgstab..rr"rr)rBTrCrDrFrFrFrF)r rGr>rIrBrJrKr r,rrrr rLr)rMr=rNr)rOrPrQrrRrSrTrVrWr+rrXrYrZr[r\r]r^r_r`rarbrcr)r=r>r?rrs`  *       z5Use Conjugate Gradient iteration to solve ``Ax = b``.zThe real or complex N-by-N matrix of the linear system. ``A`` must represent a hermitian, positive definite matrix. Alternatively, ``A`` can be a linear operator which can produce ``Ax`` using, e.g., ``scipy.sparse.linalg.LinearOperator``.c szt|||\}}}t} |dkr0| d}|j|j} tjj} tt| d} fdd} t||t j j | d}|dkr|dfS|}d}d}t d | jd }d}d}d }|}x|}| ||||||| \ }}}}}}}}|dk r||kr|t |d|d| }t |d|d| }|dkrZ|dk rV|Pn|dkr|||9<|||||7<n|d kr| ||||<n|d kr|||9<|||7<n`|d krB|rd}d}t|||\}}|dkrB|dkrB||<t|||\}}d }qW|dkrn||krn||k rn|}||fS)Nr<ZcgrevcomcstjjS)N)rrrr)r=r>r?rrr@-szcg..rr"rrr)rBTrCrDFrFrFrF)r rGr>rIrBrJrKr r,rrrr rLr)rMr=rNr)rOrPrQrrRrSrTrVrWr+rrXrYrZr[r\r]r^r_r`rarbrcr)r=r>r?rrsf  *       z=Use Conjugate Gradient Squared iteration to solve ``Ax = b``.zThe real-valued N-by-N matrix of the linear system. Alternatively, ``A`` can be a linear operator which can produce ``Ax`` using, e.g., ``scipy.sparse.linalg.LinearOperator``.c st|||\}}}t} |dkr0| d}|j|j} tjj} tt| d} fdd} t||t j j | d}|dkr|dfS|}d}d}t d | jd }d}d}d }|}x|}| ||||||| \ }}}}}}}}|dk r||kr|t |d|d| }t |d|d| }|dkrZ|dk rV|Pn|dkr|||9<|||||7<n|d kr| ||||<n|d kr|||9<|||7<n`|dkrB|rd}d}t|||\}}|dkrB|dkrB||<t|||\}}d }qW|dkrtt|\}}|rtd}|dkr||kr||k r|}||fS)Nr<Z cgsrevcomcstjjS)N)rrrr)r=r>r?rrr@tszcgs..rr"rrrd)rBTrCrDrFrFrFrFi)r rGr>rIrBrJrKr r,rrrr rLr)rMr=rNr)rOrPrQrrRrSrTrVrWr+rrXrYrZr[r\r]r^r_r`rarbrcokr)r=r>r?rrbsn  *        c ( s|dkr|}n|dk rtd|dk r>| dkr>tjdtdd| dkrJd} | dkr`td j| |dkrld } t|||\}}} t} |dkr| d }|dkrd }t|| }|j|j} t j j }t t |d }tjj}tjj| }fdd}t|| ||d} | dkr*| dfS|dkr@| dfSd}|t|| |}tj}tj}d}d}td|| j d}t|dd|dj d}d}d}d}|}|}d} d}!d}"x.|}#||||||||||| \ }}}}}}$}%}| dkr||#kr|t|d|d| }&t|d|d| }'|d kr~| d!krf|!rx|||n| dkrx|PnJ|dkr||'|%9<||'|$7<n|dkr| ||'||&<| r|dkrd}!d} n|dkrJ||'|%9<||'|$||&7<|!r| d"kr<|||d}!|"d}"n~|dkr|rbd#}d}t||&| \}}|s||krtdd|}ntdd|}|dkr|t|| |}n||}|}d}| dkr|"|kr|}PqW|dkr || k r |}| |fS)$ag Use Generalized Minimal RESidual iteration to solve ``Ax = b``. Parameters ---------- A : {sparse matrix, dense matrix, LinearOperator} The real or complex N-by-N matrix of the linear system. Alternatively, ``A`` can be a linear operator which can produce ``Ax`` using, e.g., ``scipy.sparse.linalg.LinearOperator``. 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, atol : float, optional Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``. The default for ``atol`` is ``'legacy'``, which emulates a different legacy behavior. .. warning:: The default value for `atol` will be changed in a future release. For future compatibility, specify `atol` explicitly. 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. 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(args)`, where `args` are selected by `callback_type`. callback_type : {'x', 'pr_norm', 'legacy'}, optional Callback function argument requested: - ``x``: current iterate (ndarray), called on every restart - ``pr_norm``: relative (preconditioned) residual norm (float), called on every inner iteration - ``legacy`` (default): same as ``pr_norm``, but also changes the meaning of 'maxiter' to count inner iterations instead of restart cycles. 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) Examples -------- >>> from scipy.sparse import csc_matrix >>> from scipy.sparse.linalg import gmres >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float) >>> b = np.array([2, 4, -1], dtype=float) >>> x, exitCode = gmres(A, b) >>> print(exitCode) # 0 indicates successful convergence 0 >>> np.allclose(A.dot(x), b) True NzOCannot specify both restart and restrt keywords. Preferably use 'restart' only.a4scipy.sparse.linalg.gmres called without specifying `callback_type`. The default value will be changed in a future release. For compatibility, specify a value for `callback_type` explicitly, e.g., ``{name}(..., callback_type='pr_norm')``, or to retain the old behavior ``{name}(..., callback_type='legacy')``rD)rr r!r?pr_normzUnknown callback_type: {!r}noner<Z gmresrevcomcstjjS)N)rrrr)r=r>r?rrr@:szgmres..rr"rg?rrA)rBrCTFrg?gؗҜrIrBrJrKr rrrr,nanr rLrr()(rMr=rNr)ZrestartrOrPrQZrestrtrZ callback_typerRrSrTrVrWr*ZMb_nrm2r+Zptol_max_factorZptolrZpresidrXrYrZZwork2r[r\r]r^Zold_ijobZ first_passZ resid_readyZiter_numr_r`rarbrcr)r=r>r?rrsa        0             c ! sptd|\} } |dkr|dkrtdrfdd} fdd} fdd} fd d }tj| | d }tj| |d }n(d d }tj||d }tj||d }t}|dkr|d}tjj}tt |d}fdd}t ||t j j |d}|dkr | dfS|}d}d }td|j}d}d}d}|}x|}|||||||| \ }}}}}}}}|dk r||kr|t|d|d|}t|d|d|} |d!kr|dk r|PnN|dkr$|| |9<|| |j||7<n|dkr^|| |9<|| |j||7<n|dkr||j|| ||<n|dkr|j|| ||<n|dkr|j|| ||<n~|dkr|j|| ||<n`|dkr || |9<|| |j7<n*|dkr6|r$d"}d}t|||\}}d}qPW|dkrd||krd||k rd|}| |fS)#a Use Quasi-Minimal Residual iteration to solve ``Ax = b``. Parameters ---------- A : {sparse matrix, dense matrix, LinearOperator} The real-valued N-by-N matrix of the linear system. Alternatively, ``A`` can be a linear operator which can produce ``Ax`` and ``A^T x`` using, e.g., ``scipy.sparse.linalg.LinearOperator``. 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, atol : float, optional Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``. The default for ``atol`` is ``'legacy'``, which emulates a different legacy behavior. .. warning:: The default value for `atol` will be changed in a future release. For future compatibility, specify `atol` explicitly. 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. See Also -------- LinearOperator Examples -------- >>> from scipy.sparse import csc_matrix >>> from scipy.sparse.linalg import qmr >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float) >>> b = np.array([2, 4, -1], dtype=float) >>> x, exitCode = qmr(A, b) >>> print(exitCode) # 0 indicates successful convergence 0 >>> np.allclose(A.dot(x), b) True NrTcs j|dS)Nleft)rT)r=)A_rr left_psolveszqmr..left_psolvecs j|dS)Nright)rT)r=)rmrr right_psolveszqmr..right_psolvecs j|dS)Nrl)rU)r=)rmrr left_rpsolveszqmr..left_rpsolvecs j|dS)Nro)rU)r=)rmrr right_rpsolveszqmr..right_rpsolve)r>rHcSs|S)Nr)r=rrridszqmr..idr<Z qmrrevcomcstjjjS)N)rrrr>r)rMr=r?rrr@szqmr..rr"rr TrCrDrrErArdFrFrFrF)r hasattrr shaperGrIrBrJrKr r,rrrr rLr>rHr)!rMr=rNr)rOZM1ZM2rQrrPrRrnrprqrrrsrSrVrWr+rrXrYrZr[r\r]r^r_r`rarbrcr)rMrmr=r?rrsC        *   "          )r-r.)Nr;NNNN)Nr;NNNN)Nr;NNNN)Nr;NNNN) Nr;NNNNNNN)Nr;NNNNN)r4__all__r#Znumpyrr-r Zscipy.sparse.linalg.interfacer utilsr Zscipy._lib._utilr Zscipy._lib._threadsafetyr rIr1r3rr,r:rrrrrrrrrrsL     ) ) A<AG h