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