Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_isolve/tfqmr.py: 8%
72 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-23 06:43 +0000
1import numpy as np
2from .utils import make_system
3from scipy._lib.deprecation import _deprecate_positional_args
6__all__ = ['tfqmr']
9@_deprecate_positional_args(version="1.14.0")
10def tfqmr(A, b, x0=None, *, tol=1e-5, maxiter=None, M=None,
11 callback=None, atol=None, show=False):
12 """
13 Use Transpose-Free Quasi-Minimal Residual iteration to solve ``Ax = b``.
15 Parameters
16 ----------
17 A : {sparse matrix, ndarray, LinearOperator}
18 The real or complex N-by-N matrix of the linear system.
19 Alternatively, `A` can be a linear operator which can
20 produce ``Ax`` using, e.g.,
21 `scipy.sparse.linalg.LinearOperator`.
22 b : {ndarray}
23 Right hand side of the linear system. Has shape (N,) or (N,1).
24 x0 : {ndarray}
25 Starting guess for the solution.
26 tol, atol : float, optional
27 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b-Ax0), atol)``.
28 The default for `tol` is 1.0e-5.
29 The default for `atol` is ``tol * norm(b-Ax0)``.
31 .. warning::
33 The default value for `atol` will be changed in a future release.
34 For future compatibility, specify `atol` explicitly.
35 maxiter : int, optional
36 Maximum number of iterations. Iteration will stop after maxiter
37 steps even if the specified tolerance has not been achieved.
38 Default is ``min(10000, ndofs * 10)``, where ``ndofs = A.shape[0]``.
39 M : {sparse matrix, ndarray, LinearOperator}
40 Inverse of the preconditioner of A. M should approximate the
41 inverse of A and be easy to solve for (see Notes). Effective
42 preconditioning dramatically improves the rate of convergence,
43 which implies that fewer iterations are needed to reach a given
44 error tolerance. By default, no preconditioner is used.
45 callback : function, optional
46 User-supplied function to call after each iteration. It is called
47 as `callback(xk)`, where `xk` is the current solution vector.
48 show : bool, optional
49 Specify ``show = True`` to show the convergence, ``show = False`` is
50 to close the output of the convergence.
51 Default is `False`.
53 Returns
54 -------
55 x : ndarray
56 The converged solution.
57 info : int
58 Provides convergence information:
60 - 0 : successful exit
61 - >0 : convergence to tolerance not achieved, number of iterations
62 - <0 : illegal input or breakdown
64 Notes
65 -----
66 The Transpose-Free QMR algorithm is derived from the CGS algorithm.
67 However, unlike CGS, the convergence curves for the TFQMR method is
68 smoothed by computing a quasi minimization of the residual norm. The
69 implementation supports left preconditioner, and the "residual norm"
70 to compute in convergence criterion is actually an upper bound on the
71 actual residual norm ``||b - Axk||``.
73 References
74 ----------
75 .. [1] R. W. Freund, A Transpose-Free Quasi-Minimal Residual Algorithm for
76 Non-Hermitian Linear Systems, SIAM J. Sci. Comput., 14(2), 470-482,
77 1993.
78 .. [2] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition,
79 SIAM, Philadelphia, 2003.
80 .. [3] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations,
81 number 16 in Frontiers in Applied Mathematics, SIAM, Philadelphia,
82 1995.
84 Examples
85 --------
86 >>> import numpy as np
87 >>> from scipy.sparse import csc_matrix
88 >>> from scipy.sparse.linalg import tfqmr
89 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
90 >>> b = np.array([2, 4, -1], dtype=float)
91 >>> x, exitCode = tfqmr(A, b)
92 >>> print(exitCode) # 0 indicates successful convergence
93 0
94 >>> np.allclose(A.dot(x), b)
95 True
96 """
98 # Check data type
99 dtype = A.dtype
100 if np.issubdtype(dtype, np.int64):
101 dtype = float
102 A = A.astype(dtype)
103 if np.issubdtype(b.dtype, np.int64):
104 b = b.astype(dtype)
106 A, M, x, b, postprocess = make_system(A, M, x0, b)
108 # Check if the R.H.S is a zero vector
109 if np.linalg.norm(b) == 0.:
110 x = b.copy()
111 return (postprocess(x), 0)
113 ndofs = A.shape[0]
114 if maxiter is None:
115 maxiter = min(10000, ndofs * 10)
117 if x0 is None:
118 r = b.copy()
119 else:
120 r = b - A.matvec(x)
121 u = r
122 w = r.copy()
123 # Take rstar as b - Ax0, that is rstar := r = b - Ax0 mathematically
124 rstar = r
125 v = M.matvec(A.matvec(r))
126 uhat = v
127 d = theta = eta = 0.
128 rho = np.inner(rstar.conjugate(), r)
129 rhoLast = rho
130 r0norm = np.sqrt(rho)
131 tau = r0norm
132 if r0norm == 0:
133 return (postprocess(x), 0)
135 if atol is None:
136 atol = tol * r0norm
137 else:
138 atol = max(atol, tol * r0norm)
140 for iter in range(maxiter):
141 even = iter % 2 == 0
142 if (even):
143 vtrstar = np.inner(rstar.conjugate(), v)
144 # Check breakdown
145 if vtrstar == 0.:
146 return (postprocess(x), -1)
147 alpha = rho / vtrstar
148 uNext = u - alpha * v # [1]-(5.6)
149 w -= alpha * uhat # [1]-(5.8)
150 d = u + (theta**2 / alpha) * eta * d # [1]-(5.5)
151 # [1]-(5.2)
152 theta = np.linalg.norm(w) / tau
153 c = np.sqrt(1. / (1 + theta**2))
154 tau *= theta * c
155 # Calculate step and direction [1]-(5.4)
156 eta = (c**2) * alpha
157 z = M.matvec(d)
158 x += eta * z
160 if callback is not None:
161 callback(x)
163 # Convergence criteron
164 if tau * np.sqrt(iter+1) < atol:
165 if (show):
166 print("TFQMR: Linear solve converged due to reach TOL "
167 "iterations {}".format(iter+1))
168 return (postprocess(x), 0)
170 if (not even):
171 # [1]-(5.7)
172 rho = np.inner(rstar.conjugate(), w)
173 beta = rho / rhoLast
174 u = w + beta * u
175 v = beta * uhat + (beta**2) * v
176 uhat = M.matvec(A.matvec(u))
177 v += uhat
178 else:
179 uhat = M.matvec(A.matvec(uNext))
180 u = uNext
181 rhoLast = rho
183 if (show):
184 print("TFQMR: Linear solve not converged due to reach MAXIT "
185 "iterations {}".format(iter+1))
186 return (postprocess(x), maxiter)