import numpy as np
from .iterative import _get_atol_rtol
from .utils import make_system
__all__ = ['tfqmr']
[文档]
def tfqmr(A, b, x0=None, *, rtol=1e-5, atol=0., maxiter=None, M=None,
callback=None, show=False):
"""
Use Transpose-Free Quasi-Minimal Residual iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse matrix, ndarray, 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 : {ndarray}
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : {ndarray}
Starting guess for the solution.
rtol, atol : float, optional
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``rtol=1e-5``, the default for ``atol`` is ``0.0``.
maxiter : int, optional
Maximum number of iterations. Iteration will stop after maxiter
steps even if the specified tolerance has not been achieved.
Default is ``min(10000, ndofs * 10)``, where ``ndofs = A.shape[0]``.
M : {sparse matrix, ndarray, 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, optional
User-supplied function to call after each iteration. It is called
as `callback(xk)`, where `xk` is the current solution vector.
show : bool, optional
Specify ``show = True`` to show the convergence, ``show = False`` is
to close the output of the convergence.
Default is `False`.
Returns
-------
x : ndarray
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
Notes
-----
The Transpose-Free QMR algorithm is derived from the CGS algorithm.
However, unlike CGS, the convergence curves for the TFQMR method is
smoothed by computing a quasi minimization of the residual norm. The
implementation supports left preconditioner, and the "residual norm"
to compute in convergence criterion is actually an upper bound on the
actual residual norm ``||b - Axk||``.
References
----------
.. [1] R. W. Freund, A Transpose-Free Quasi-Minimal Residual Algorithm for
Non-Hermitian Linear Systems, SIAM J. Sci. Comput., 14(2), 470-482,
1993.
.. [2] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition,
SIAM, Philadelphia, 2003.
.. [3] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations,
number 16 in Frontiers in Applied Mathematics, SIAM, Philadelphia,
1995.
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import tfqmr
>>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
>>> b = np.array([2, 4, -1], dtype=float)
>>> x, exitCode = tfqmr(A, b, atol=0.0)
>>> print(exitCode) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
# Check data type
dtype = A.dtype
if np.issubdtype(dtype, np.int64):
dtype = float
A = A.astype(dtype)
if np.issubdtype(b.dtype, np.int64):
b = b.astype(dtype)
A, M, x, b, postprocess = make_system(A, M, x0, b)
# Check if the R.H.S is a zero vector
if np.linalg.norm(b) == 0.:
x = b.copy()
return (postprocess(x), 0)
ndofs = A.shape[0]
if maxiter is None:
maxiter = min(10000, ndofs * 10)
if x0 is None:
r = b.copy()
else:
r = b - A.matvec(x)
u = r
w = r.copy()
# Take rstar as b - Ax0, that is rstar := r = b - Ax0 mathematically
rstar = r
v = M.matvec(A.matvec(r))
uhat = v
d = theta = eta = 0.
# at this point we know rstar == r, so rho is always real
rho = np.inner(rstar.conjugate(), r).real
rhoLast = rho
r0norm = np.sqrt(rho)
tau = r0norm
if r0norm == 0:
return (postprocess(x), 0)
# we call this to get the right atol and raise errors as necessary
atol, _ = _get_atol_rtol('tfqmr', r0norm, atol, rtol)
for iter in range(maxiter):
even = iter % 2 == 0
if (even):
vtrstar = np.inner(rstar.conjugate(), v)
# Check breakdown
if vtrstar == 0.:
return (postprocess(x), -1)
alpha = rho / vtrstar
uNext = u - alpha * v # [1]-(5.6)
w -= alpha * uhat # [1]-(5.8)
d = u + (theta**2 / alpha) * eta * d # [1]-(5.5)
# [1]-(5.2)
theta = np.linalg.norm(w) / tau
c = np.sqrt(1. / (1 + theta**2))
tau *= theta * c
# Calculate step and direction [1]-(5.4)
eta = (c**2) * alpha
z = M.matvec(d)
x += eta * z
if callback is not None:
callback(x)
# Convergence criterion
if tau * np.sqrt(iter+1) < atol:
if (show):
print("TFQMR: Linear solve converged due to reach TOL "
f"iterations {iter+1}")
return (postprocess(x), 0)
if (not even):
# [1]-(5.7)
rho = np.inner(rstar.conjugate(), w)
beta = rho / rhoLast
u = w + beta * u
v = beta * uhat + (beta**2) * v
uhat = M.matvec(A.matvec(u))
v += uhat
else:
uhat = M.matvec(A.matvec(uNext))
u = uNext
rhoLast = rho
if (show):
print("TFQMR: Linear solve not converged due to reach MAXIT "
f"iterations {iter+1}")
return (postprocess(x), maxiter)
