Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_isolve/minres.py: 4%
203 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-03-22 06:44 +0000
1import warnings
2from numpy import inner, zeros, inf, finfo
3from numpy.linalg import norm
4from math import sqrt
6from .utils import make_system
7from scipy._lib.deprecation import _NoValue, _deprecate_positional_args
9__all__ = ['minres']
12@_deprecate_positional_args(version="1.14.0")
13def minres(A, b, x0=None, *, shift=0.0, tol=_NoValue, maxiter=None,
14 M=None, callback=None, show=False, check=False, rtol=1e-5):
15 """
16 Use MINimum RESidual iteration to solve Ax=b
18 MINRES minimizes norm(Ax - b) for a real symmetric matrix A. Unlike
19 the Conjugate Gradient method, A can be indefinite or singular.
21 If shift != 0 then the method solves (A - shift*I)x = b
23 Parameters
24 ----------
25 A : {sparse matrix, ndarray, LinearOperator}
26 The real symmetric N-by-N matrix of the linear system
27 Alternatively, ``A`` can be a linear operator which can
28 produce ``Ax`` using, e.g.,
29 ``scipy.sparse.linalg.LinearOperator``.
30 b : ndarray
31 Right hand side of the linear system. Has shape (N,) or (N,1).
33 Returns
34 -------
35 x : ndarray
36 The converged solution.
37 info : integer
38 Provides convergence information:
39 0 : successful exit
40 >0 : convergence to tolerance not achieved, number of iterations
41 <0 : illegal input or breakdown
43 Other Parameters
44 ----------------
45 x0 : ndarray
46 Starting guess for the solution.
47 shift : float
48 Value to apply to the system ``(A - shift * I)x = b``. Default is 0.
49 rtol : float
50 Tolerance to achieve. The algorithm terminates when the relative
51 residual is below ``rtol``.
52 maxiter : integer
53 Maximum number of iterations. Iteration will stop after maxiter
54 steps even if the specified tolerance has not been achieved.
55 M : {sparse matrix, ndarray, LinearOperator}
56 Preconditioner for A. The preconditioner should approximate the
57 inverse of A. Effective preconditioning dramatically improves the
58 rate of convergence, which implies that fewer iterations are needed
59 to reach a given error tolerance.
60 callback : function
61 User-supplied function to call after each iteration. It is called
62 as callback(xk), where xk is the current solution vector.
63 show : bool
64 If ``True``, print out a summary and metrics related to the solution
65 during iterations. Default is ``False``.
66 check : bool
67 If ``True``, run additional input validation to check that `A` and
68 `M` (if specified) are symmetric. Default is ``False``.
69 tol : float, optional, deprecated
71 .. deprecated:: 1.12.0
72 `minres` keyword argument ``tol`` is deprecated in favor of ``rtol``
73 and will be removed in SciPy 1.14.0.
75 Examples
76 --------
77 >>> import numpy as np
78 >>> from scipy.sparse import csc_matrix
79 >>> from scipy.sparse.linalg import minres
80 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
81 >>> A = A + A.T
82 >>> b = np.array([2, 4, -1], dtype=float)
83 >>> x, exitCode = minres(A, b)
84 >>> print(exitCode) # 0 indicates successful convergence
85 0
86 >>> np.allclose(A.dot(x), b)
87 True
89 References
90 ----------
91 Solution of sparse indefinite systems of linear equations,
92 C. C. Paige and M. A. Saunders (1975),
93 SIAM J. Numer. Anal. 12(4), pp. 617-629.
94 https://web.stanford.edu/group/SOL/software/minres/
96 This file is a translation of the following MATLAB implementation:
97 https://web.stanford.edu/group/SOL/software/minres/minres-matlab.zip
99 """
100 A, M, x, b, postprocess = make_system(A, M, x0, b)
102 if tol is not _NoValue:
103 msg = ("'scipy.sparse.linalg.minres' keyword argument `tol` is "
104 "deprecated in favor of `rtol` and will be removed in SciPy "
105 "v1.14. Until then, if set, it will override `rtol`.")
106 warnings.warn(msg, category=DeprecationWarning, stacklevel=4)
107 rtol = float(tol) if tol is not None else rtol
109 matvec = A.matvec
110 psolve = M.matvec
112 first = 'Enter minres. '
113 last = 'Exit minres. '
115 n = A.shape[0]
117 if maxiter is None:
118 maxiter = 5 * n
120 msg = [' beta2 = 0. If M = I, b and x are eigenvectors ', # -1
121 ' beta1 = 0. The exact solution is x0 ', # 0
122 ' A solution to Ax = b was found, given rtol ', # 1
123 ' A least-squares solution was found, given rtol ', # 2
124 ' Reasonable accuracy achieved, given eps ', # 3
125 ' x has converged to an eigenvector ', # 4
126 ' acond has exceeded 0.1/eps ', # 5
127 ' The iteration limit was reached ', # 6
128 ' A does not define a symmetric matrix ', # 7
129 ' M does not define a symmetric matrix ', # 8
130 ' M does not define a pos-def preconditioner '] # 9
132 if show:
133 print(first + 'Solution of symmetric Ax = b')
134 print(first + f'n = {n:3g} shift = {shift:23.14e}')
135 print(first + f'itnlim = {maxiter:3g} rtol = {rtol:11.2e}')
136 print()
138 istop = 0
139 itn = 0
140 Anorm = 0
141 Acond = 0
142 rnorm = 0
143 ynorm = 0
145 xtype = x.dtype
147 eps = finfo(xtype).eps
149 # Set up y and v for the first Lanczos vector v1.
150 # y = beta1 P' v1, where P = C**(-1).
151 # v is really P' v1.
153 if x0 is None:
154 r1 = b.copy()
155 else:
156 r1 = b - A@x
157 y = psolve(r1)
159 beta1 = inner(r1, y)
161 if beta1 < 0:
162 raise ValueError('indefinite preconditioner')
163 elif beta1 == 0:
164 return (postprocess(x), 0)
166 bnorm = norm(b)
167 if bnorm == 0:
168 x = b
169 return (postprocess(x), 0)
171 beta1 = sqrt(beta1)
173 if check:
174 # are these too strict?
176 # see if A is symmetric
177 w = matvec(y)
178 r2 = matvec(w)
179 s = inner(w,w)
180 t = inner(y,r2)
181 z = abs(s - t)
182 epsa = (s + eps) * eps**(1.0/3.0)
183 if z > epsa:
184 raise ValueError('non-symmetric matrix')
186 # see if M is symmetric
187 r2 = psolve(y)
188 s = inner(y,y)
189 t = inner(r1,r2)
190 z = abs(s - t)
191 epsa = (s + eps) * eps**(1.0/3.0)
192 if z > epsa:
193 raise ValueError('non-symmetric preconditioner')
195 # Initialize other quantities
196 oldb = 0
197 beta = beta1
198 dbar = 0
199 epsln = 0
200 qrnorm = beta1
201 phibar = beta1
202 rhs1 = beta1
203 rhs2 = 0
204 tnorm2 = 0
205 gmax = 0
206 gmin = finfo(xtype).max
207 cs = -1
208 sn = 0
209 w = zeros(n, dtype=xtype)
210 w2 = zeros(n, dtype=xtype)
211 r2 = r1
213 if show:
214 print()
215 print()
216 print(' Itn x(1) Compatible LS norm(A) cond(A) gbar/|A|')
218 while itn < maxiter:
219 itn += 1
221 s = 1.0/beta
222 v = s*y
224 y = matvec(v)
225 y = y - shift * v
227 if itn >= 2:
228 y = y - (beta/oldb)*r1
230 alfa = inner(v,y)
231 y = y - (alfa/beta)*r2
232 r1 = r2
233 r2 = y
234 y = psolve(r2)
235 oldb = beta
236 beta = inner(r2,y)
237 if beta < 0:
238 raise ValueError('non-symmetric matrix')
239 beta = sqrt(beta)
240 tnorm2 += alfa**2 + oldb**2 + beta**2
242 if itn == 1:
243 if beta/beta1 <= 10*eps:
244 istop = -1 # Terminate later
246 # Apply previous rotation Qk-1 to get
247 # [deltak epslnk+1] = [cs sn][dbark 0 ]
248 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
250 oldeps = epsln
251 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
252 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
253 epsln = sn * beta # epsln2 = 0 epslnk+1
254 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
255 root = norm([gbar, dbar])
256 Arnorm = phibar * root
258 # Compute the next plane rotation Qk
260 gamma = norm([gbar, beta]) # gammak
261 gamma = max(gamma, eps)
262 cs = gbar / gamma # ck
263 sn = beta / gamma # sk
264 phi = cs * phibar # phik
265 phibar = sn * phibar # phibark+1
267 # Update x.
269 denom = 1.0/gamma
270 w1 = w2
271 w2 = w
272 w = (v - oldeps*w1 - delta*w2) * denom
273 x = x + phi*w
275 # Go round again.
277 gmax = max(gmax, gamma)
278 gmin = min(gmin, gamma)
279 z = rhs1 / gamma
280 rhs1 = rhs2 - delta*z
281 rhs2 = - epsln*z
283 # Estimate various norms and test for convergence.
285 Anorm = sqrt(tnorm2)
286 ynorm = norm(x)
287 epsa = Anorm * eps
288 epsx = Anorm * ynorm * eps
289 epsr = Anorm * ynorm * rtol
290 diag = gbar
292 if diag == 0:
293 diag = epsa
295 qrnorm = phibar
296 rnorm = qrnorm
297 if ynorm == 0 or Anorm == 0:
298 test1 = inf
299 else:
300 test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||)
301 if Anorm == 0:
302 test2 = inf
303 else:
304 test2 = root / Anorm # ||Ar|| / (||A|| ||r||)
306 # Estimate cond(A).
307 # In this version we look at the diagonals of R in the
308 # factorization of the lower Hessenberg matrix, Q @ H = R,
309 # where H is the tridiagonal matrix from Lanczos with one
310 # extra row, beta(k+1) e_k^T.
312 Acond = gmax/gmin
314 # See if any of the stopping criteria are satisfied.
315 # In rare cases, istop is already -1 from above (Abar = const*I).
317 if istop == 0:
318 t1 = 1 + test1 # These tests work if rtol < eps
319 t2 = 1 + test2
320 if t2 <= 1:
321 istop = 2
322 if t1 <= 1:
323 istop = 1
325 if itn >= maxiter:
326 istop = 6
327 if Acond >= 0.1/eps:
328 istop = 4
329 if epsx >= beta1:
330 istop = 3
331 # if rnorm <= epsx : istop = 2
332 # if rnorm <= epsr : istop = 1
333 if test2 <= rtol:
334 istop = 2
335 if test1 <= rtol:
336 istop = 1
338 # See if it is time to print something.
340 prnt = False
341 if n <= 40:
342 prnt = True
343 if itn <= 10:
344 prnt = True
345 if itn >= maxiter-10:
346 prnt = True
347 if itn % 10 == 0:
348 prnt = True
349 if qrnorm <= 10*epsx:
350 prnt = True
351 if qrnorm <= 10*epsr:
352 prnt = True
353 if Acond <= 1e-2/eps:
354 prnt = True
355 if istop != 0:
356 prnt = True
358 if show and prnt:
359 str1 = f'{itn:6g} {x[0]:12.5e} {test1:10.3e}'
360 str2 = f' {test2:10.3e}'
361 str3 = f' {Anorm:8.1e} {Acond:8.1e} {gbar/Anorm:8.1e}'
363 print(str1 + str2 + str3)
365 if itn % 10 == 0:
366 print()
368 if callback is not None:
369 callback(x)
371 if istop != 0:
372 break # TODO check this
374 if show:
375 print()
376 print(last + f' istop = {istop:3g} itn ={itn:5g}')
377 print(last + f' Anorm = {Anorm:12.4e} Acond = {Acond:12.4e}')
378 print(last + f' rnorm = {rnorm:12.4e} ynorm = {ynorm:12.4e}')
379 print(last + f' Arnorm = {Arnorm:12.4e}')
380 print(last + msg[istop+1])
382 if istop == 6:
383 info = maxiter
384 else:
385 info = 0
387 return (postprocess(x),info)