Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_isolve/tfqmr.py: 10%

1import numpy as np 

2from .iterative import _get_atol_rtol 

3from .utils import make_system 

4from scipy._lib.deprecation import _NoValue, _deprecate_positional_args 

5 

6 

7__all__ = ['tfqmr'] 

8 

9 

10@_deprecate_positional_args(version="1.14.0") 

11def tfqmr(A, b, x0=None, *, tol=_NoValue, maxiter=None, M=None, 

12 callback=None, atol=None, rtol=1e-5, show=False): 

13 """ 

14 Use Transpose-Free Quasi-Minimal Residual iteration to solve ``Ax = b``. 

15 

16 Parameters 

17 ---------- 

18 A : {sparse matrix, ndarray, LinearOperator} 

19 The real or complex N-by-N matrix of the linear system. 

20 Alternatively, `A` can be a linear operator which can 

21 produce ``Ax`` using, e.g., 

22 `scipy.sparse.linalg.LinearOperator`. 

23 b : {ndarray} 

24 Right hand side of the linear system. Has shape (N,) or (N,1). 

25 x0 : {ndarray} 

26 Starting guess for the solution. 

27 rtol, atol : float, optional 

28 Parameters for the convergence test. For convergence, 

29 ``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied. 

30 The default is ``rtol=1e-5``, the default for ``atol`` is ``rtol``. 

31 

32 .. warning:: 

33 

34 The default value for ``atol`` will be changed to ``0.0`` in 

35 SciPy 1.14.0. 

36 maxiter : int, optional 

37 Maximum number of iterations. Iteration will stop after maxiter 

38 steps even if the specified tolerance has not been achieved. 

39 Default is ``min(10000, ndofs * 10)``, where ``ndofs = A.shape[0]``. 

40 M : {sparse matrix, ndarray, LinearOperator} 

41 Inverse of the preconditioner of A. M should approximate the 

42 inverse of A and be easy to solve for (see Notes). Effective 

43 preconditioning dramatically improves the rate of convergence, 

44 which implies that fewer iterations are needed to reach a given 

45 error tolerance. By default, no preconditioner is used. 

46 callback : function, optional 

47 User-supplied function to call after each iteration. It is called 

48 as `callback(xk)`, where `xk` is the current solution vector. 

49 show : bool, optional 

50 Specify ``show = True`` to show the convergence, ``show = False`` is 

51 to close the output of the convergence. 

52 Default is `False`. 

53 tol : float, optional, deprecated 

54 

55 .. deprecated:: 1.12.0 

56 `tfqmr` keyword argument ``tol`` is deprecated in favor of ``rtol`` 

57 and will be removed in SciPy 1.14.0. 

58 

59 Returns 

60 ------- 

61 x : ndarray 

62 The converged solution. 

63 info : int 

64 Provides convergence information: 

65 

66 - 0 : successful exit 

67 - >0 : convergence to tolerance not achieved, number of iterations 

68 - <0 : illegal input or breakdown 

69 

70 Notes 

71 ----- 

72 The Transpose-Free QMR algorithm is derived from the CGS algorithm. 

73 However, unlike CGS, the convergence curves for the TFQMR method is 

74 smoothed by computing a quasi minimization of the residual norm. The 

75 implementation supports left preconditioner, and the "residual norm" 

76 to compute in convergence criterion is actually an upper bound on the 

77 actual residual norm ``||b - Axk||``. 

78 

79 References 

80 ---------- 

81 .. [1] R. W. Freund, A Transpose-Free Quasi-Minimal Residual Algorithm for 

82 Non-Hermitian Linear Systems, SIAM J. Sci. Comput., 14(2), 470-482, 

83 1993. 

84 .. [2] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, 

85 SIAM, Philadelphia, 2003. 

86 .. [3] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, 

87 number 16 in Frontiers in Applied Mathematics, SIAM, Philadelphia, 

88 1995. 

89 

90 Examples 

91 -------- 

92 >>> import numpy as np 

93 >>> from scipy.sparse import csc_matrix 

94 >>> from scipy.sparse.linalg import tfqmr 

95 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float) 

96 >>> b = np.array([2, 4, -1], dtype=float) 

97 >>> x, exitCode = tfqmr(A, b) 

98 >>> print(exitCode) # 0 indicates successful convergence 

99 0 

100 >>> np.allclose(A.dot(x), b) 

101 True 

102 """ 

103 

104 # Check data type 

105 dtype = A.dtype 

106 if np.issubdtype(dtype, np.int64): 

107 dtype = float 

108 A = A.astype(dtype) 

109 if np.issubdtype(b.dtype, np.int64): 

110 b = b.astype(dtype) 

111 

112 A, M, x, b, postprocess = make_system(A, M, x0, b) 

113 

114 # Check if the R.H.S is a zero vector 

115 if np.linalg.norm(b) == 0.: 

116 x = b.copy() 

117 return (postprocess(x), 0) 

118 

119 ndofs = A.shape[0] 

120 if maxiter is None: 

121 maxiter = min(10000, ndofs * 10) 

122 

123 if x0 is None: 

124 r = b.copy() 

125 else: 

126 r = b - A.matvec(x) 

127 u = r 

128 w = r.copy() 

129 # Take rstar as b - Ax0, that is rstar := r = b - Ax0 mathematically 

130 rstar = r 

131 v = M.matvec(A.matvec(r)) 

132 uhat = v 

133 d = theta = eta = 0. 

134 # at this point we know rstar == r, so rho is always real 

135 rho = np.inner(rstar.conjugate(), r).real 

136 rhoLast = rho 

137 r0norm = np.sqrt(rho) 

138 tau = r0norm 

139 if r0norm == 0: 

140 return (postprocess(x), 0) 

141 

142 # we call this to get the right atol and raise warnings as necessary 

143 atol, _ = _get_atol_rtol('tfqmr', r0norm, tol, atol, rtol) 

144 

145 for iter in range(maxiter): 

146 even = iter % 2 == 0 

147 if (even): 

148 vtrstar = np.inner(rstar.conjugate(), v) 

149 # Check breakdown 

150 if vtrstar == 0.: 

151 return (postprocess(x), -1) 

152 alpha = rho / vtrstar 

153 uNext = u - alpha * v # [1]-(5.6) 

154 w -= alpha * uhat # [1]-(5.8) 

155 d = u + (theta**2 / alpha) * eta * d # [1]-(5.5) 

156 # [1]-(5.2) 

157 theta = np.linalg.norm(w) / tau 

158 c = np.sqrt(1. / (1 + theta**2)) 

159 tau *= theta * c 

160 # Calculate step and direction [1]-(5.4) 

161 eta = (c**2) * alpha 

162 z = M.matvec(d) 

163 x += eta * z 

164 

165 if callback is not None: 

166 callback(x) 

167 

168 # Convergence criterion 

169 if tau * np.sqrt(iter+1) < atol: 

170 if (show): 

171 print("TFQMR: Linear solve converged due to reach TOL " 

172 f"iterations {iter+1}") 

173 return (postprocess(x), 0) 

174 

175 if (not even): 

176 # [1]-(5.7) 

177 rho = np.inner(rstar.conjugate(), w) 

178 beta = rho / rhoLast 

179 u = w + beta * u 

180 v = beta * uhat + (beta**2) * v 

181 uhat = M.matvec(A.matvec(u)) 

182 v += uhat 

183 else: 

184 uhat = M.matvec(A.matvec(uNext)) 

185 u = uNext 

186 rhoLast = rho 

187 

188 if (show): 

189 print("TFQMR: Linear solve not converged due to reach MAXIT " 

190 f"iterations {iter+1}") 

191 return (postprocess(x), maxiter)