Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/iterative.py: 8%
423 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Iterative methods for solving linear systems"""
3__all__ = ['bicg','bicgstab','cg','cgs','gmres','qmr']
5import warnings
6from textwrap import dedent
7import numpy as np
9from . import _iterative
11from scipy.sparse.linalg._interface import LinearOperator
12from .utils import make_system
13from scipy._lib._util import _aligned_zeros
14from scipy._lib._threadsafety import non_reentrant
16_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
19# Part of the docstring common to all iterative solvers
20common_doc1 = \
21"""
22Parameters
23----------
24A : {sparse matrix, ndarray, LinearOperator}"""
26common_doc2 = \
27"""b : ndarray
28 Right hand side of the linear system. Has shape (N,) or (N,1).
30Returns
31-------
32x : ndarray
33 The converged solution.
34info : integer
35 Provides convergence information:
36 0 : successful exit
37 >0 : convergence to tolerance not achieved, number of iterations
38 <0 : illegal input or breakdown
40Other Parameters
41----------------
42x0 : ndarray
43 Starting guess for the solution.
44tol, atol : float, optional
45 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``.
46 The default for ``atol`` is ``'legacy'``, which emulates
47 a different legacy behavior.
49 .. warning::
51 The default value for `atol` will be changed in a future release.
52 For future compatibility, specify `atol` explicitly.
53maxiter : integer
54 Maximum number of iterations. Iteration will stop after maxiter
55 steps even if the specified tolerance has not been achieved.
56M : {sparse matrix, ndarray, LinearOperator}
57 Preconditioner for A. The preconditioner should approximate the
58 inverse of A. Effective preconditioning dramatically improves the
59 rate of convergence, which implies that fewer iterations are needed
60 to reach a given error tolerance.
61callback : function
62 User-supplied function to call after each iteration. It is called
63 as callback(xk), where xk is the current solution vector.
64"""
67def _stoptest(residual, atol):
68 """
69 Successful termination condition for the solvers.
70 """
71 resid = np.linalg.norm(residual)
72 if resid <= atol:
73 return resid, 1
74 else:
75 return resid, 0
78def _get_atol(tol, atol, bnrm2, get_residual, routine_name):
79 """
80 Parse arguments for absolute tolerance in termination condition.
82 Parameters
83 ----------
84 tol, atol : object
85 The arguments passed into the solver routine by user.
86 bnrm2 : float
87 2-norm of the rhs vector.
88 get_residual : callable
89 Callable ``get_residual()`` that returns the initial value of
90 the residual.
91 routine_name : str
92 Name of the routine.
93 """
95 if atol is None:
96 warnings.warn("scipy.sparse.linalg.{name} called without specifying `atol`. "
97 "The default value will be changed in a future release. "
98 "For compatibility, specify a value for `atol` explicitly, e.g., "
99 "``{name}(..., atol=0)``, or to retain the old behavior "
100 "``{name}(..., atol='legacy')``".format(name=routine_name),
101 category=DeprecationWarning, stacklevel=4)
102 atol = 'legacy'
104 tol = float(tol)
106 if atol == 'legacy':
107 # emulate old legacy behavior
108 resid = get_residual()
109 if resid <= tol:
110 return 'exit'
111 if bnrm2 == 0:
112 return tol
113 else:
114 return tol * float(bnrm2)
115 else:
116 return max(float(atol), tol * float(bnrm2))
119def set_docstring(header, Ainfo, footer='', atol_default='0'):
120 def combine(fn):
121 fn.__doc__ = '\n'.join((header, common_doc1,
122 ' ' + Ainfo.replace('\n', '\n '),
123 common_doc2, dedent(footer)))
124 return fn
125 return combine
128@set_docstring('Use BIConjugate Gradient iteration to solve ``Ax = b``.',
129 'The real or complex N-by-N matrix of the linear system.\n'
130 'Alternatively, ``A`` can be a linear operator which can\n'
131 'produce ``Ax`` and ``A^T x`` using, e.g.,\n'
132 '``scipy.sparse.linalg.LinearOperator``.',
133 footer="""\
134 Examples
135 --------
136 >>> import numpy as np
137 >>> from scipy.sparse import csc_matrix
138 >>> from scipy.sparse.linalg import bicg
139 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
140 >>> b = np.array([2, 4, -1], dtype=float)
141 >>> x, exitCode = bicg(A, b)
142 >>> print(exitCode) # 0 indicates successful convergence
143 0
144 >>> np.allclose(A.dot(x), b)
145 True
147 """
148 )
149@non_reentrant()
150def bicg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None):
151 A,M,x,b,postprocess = make_system(A, M, x0, b)
153 n = len(b)
154 if maxiter is None:
155 maxiter = n*10
157 matvec, rmatvec = A.matvec, A.rmatvec
158 psolve, rpsolve = M.matvec, M.rmatvec
159 ltr = _type_conv[x.dtype.char]
160 revcom = getattr(_iterative, ltr + 'bicgrevcom')
162 get_residual = lambda: np.linalg.norm(matvec(x) - b)
163 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'bicg')
164 if atol == 'exit':
165 return postprocess(x), 0
167 resid = atol
168 ndx1 = 1
169 ndx2 = -1
170 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
171 work = _aligned_zeros(6*n,dtype=x.dtype)
172 ijob = 1
173 info = 0
174 ftflag = True
175 iter_ = maxiter
176 while True:
177 olditer = iter_
178 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
179 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
180 if callback is not None and iter_ > olditer:
181 callback(x)
182 slice1 = slice(ndx1-1, ndx1-1+n)
183 slice2 = slice(ndx2-1, ndx2-1+n)
184 if (ijob == -1):
185 if callback is not None:
186 callback(x)
187 break
188 elif (ijob == 1):
189 work[slice2] *= sclr2
190 work[slice2] += sclr1*matvec(work[slice1])
191 elif (ijob == 2):
192 work[slice2] *= sclr2
193 work[slice2] += sclr1*rmatvec(work[slice1])
194 elif (ijob == 3):
195 work[slice1] = psolve(work[slice2])
196 elif (ijob == 4):
197 work[slice1] = rpsolve(work[slice2])
198 elif (ijob == 5):
199 work[slice2] *= sclr2
200 work[slice2] += sclr1*matvec(x)
201 elif (ijob == 6):
202 if ftflag:
203 info = -1
204 ftflag = False
205 resid, info = _stoptest(work[slice1], atol)
206 ijob = 2
208 if info > 0 and iter_ == maxiter and not (resid <= atol):
209 # info isn't set appropriately otherwise
210 info = iter_
212 return postprocess(x), info
215@set_docstring('Use BIConjugate Gradient STABilized iteration to solve '
216 '``Ax = b``.',
217 'The real or complex N-by-N matrix of the linear system.\n'
218 'Alternatively, ``A`` can be a linear operator which can\n'
219 'produce ``Ax`` using, e.g.,\n'
220 '``scipy.sparse.linalg.LinearOperator``.',
221 footer="""\
222 Examples
223 --------
224 >>> import numpy as np
225 >>> from scipy.sparse import csc_matrix
226 >>> from scipy.sparse.linalg import bicgstab
227 >>> R = np.array([[4, 2, 0, 1],
228 ... [3, 0, 0, 2],
229 ... [0, 1, 1, 1],
230 ... [0, 2, 1, 0]])
231 >>> A = csc_matrix(R)
232 >>> b = np.array([-1, -0.5, -1, 2])
233 >>> x, exit_code = bicgstab(A, b)
234 >>> print(exit_code) # 0 indicates successful convergence
235 0
236 >>> np.allclose(A.dot(x), b)
237 True
238 """)
239@non_reentrant()
240def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None):
241 A, M, x, b, postprocess = make_system(A, M, x0, b)
243 n = len(b)
244 if maxiter is None:
245 maxiter = n*10
247 matvec = A.matvec
248 psolve = M.matvec
249 ltr = _type_conv[x.dtype.char]
250 revcom = getattr(_iterative, ltr + 'bicgstabrevcom')
252 get_residual = lambda: np.linalg.norm(matvec(x) - b)
253 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'bicgstab')
254 if atol == 'exit':
255 return postprocess(x), 0
257 resid = atol
258 ndx1 = 1
259 ndx2 = -1
260 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
261 work = _aligned_zeros(7*n,dtype=x.dtype)
262 ijob = 1
263 info = 0
264 ftflag = True
265 iter_ = maxiter
266 while True:
267 olditer = iter_
268 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
269 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
270 if callback is not None and iter_ > olditer:
271 callback(x)
272 slice1 = slice(ndx1-1, ndx1-1+n)
273 slice2 = slice(ndx2-1, ndx2-1+n)
274 if (ijob == -1):
275 if callback is not None:
276 callback(x)
277 break
278 elif (ijob == 1):
279 work[slice2] *= sclr2
280 work[slice2] += sclr1*matvec(work[slice1])
281 elif (ijob == 2):
282 work[slice1] = psolve(work[slice2])
283 elif (ijob == 3):
284 work[slice2] *= sclr2
285 work[slice2] += sclr1*matvec(x)
286 elif (ijob == 4):
287 if ftflag:
288 info = -1
289 ftflag = False
290 resid, info = _stoptest(work[slice1], atol)
291 ijob = 2
293 if info > 0 and iter_ == maxiter and not (resid <= atol):
294 # info isn't set appropriately otherwise
295 info = iter_
297 return postprocess(x), info
300@set_docstring('Use Conjugate Gradient iteration to solve ``Ax = b``.',
301 'The real or complex N-by-N matrix of the linear system.\n'
302 '``A`` must represent a hermitian, positive definite matrix.\n'
303 'Alternatively, ``A`` can be a linear operator which can\n'
304 'produce ``Ax`` using, e.g.,\n'
305 '``scipy.sparse.linalg.LinearOperator``.',
306 footer="""\
307 Examples
308 --------
309 >>> import numpy as np
310 >>> from scipy.sparse import csc_matrix
311 >>> from scipy.sparse.linalg import cg
312 >>> P = np.array([[4, 0, 1, 0],
313 ... [0, 5, 0, 0],
314 ... [1, 0, 3, 2],
315 ... [0, 0, 2, 4]])
316 >>> A = csc_matrix(P)
317 >>> b = np.array([-1, -0.5, -1, 2])
318 >>> x, exit_code = cg(A, b)
319 >>> print(exit_code) # 0 indicates successful convergence
320 0
321 >>> np.allclose(A.dot(x), b)
322 True
324 """)
325@non_reentrant()
326def cg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None):
327 A, M, x, b, postprocess = make_system(A, M, x0, b)
329 n = len(b)
330 if maxiter is None:
331 maxiter = n*10
333 matvec = A.matvec
334 psolve = M.matvec
335 ltr = _type_conv[x.dtype.char]
336 revcom = getattr(_iterative, ltr + 'cgrevcom')
338 get_residual = lambda: np.linalg.norm(matvec(x) - b)
339 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'cg')
340 if atol == 'exit':
341 return postprocess(x), 0
343 resid = atol
344 ndx1 = 1
345 ndx2 = -1
346 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
347 work = _aligned_zeros(4*n,dtype=x.dtype)
348 ijob = 1
349 info = 0
350 ftflag = True
351 iter_ = maxiter
352 while True:
353 olditer = iter_
354 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
355 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
356 if callback is not None and iter_ > olditer:
357 callback(x)
358 slice1 = slice(ndx1-1, ndx1-1+n)
359 slice2 = slice(ndx2-1, ndx2-1+n)
360 if (ijob == -1):
361 if callback is not None:
362 callback(x)
363 break
364 elif (ijob == 1):
365 work[slice2] *= sclr2
366 work[slice2] += sclr1*matvec(work[slice1])
367 elif (ijob == 2):
368 work[slice1] = psolve(work[slice2])
369 elif (ijob == 3):
370 work[slice2] *= sclr2
371 work[slice2] += sclr1*matvec(x)
372 elif (ijob == 4):
373 if ftflag:
374 info = -1
375 ftflag = False
376 resid, info = _stoptest(work[slice1], atol)
377 if info == 1 and iter_ > 1:
378 # recompute residual and recheck, to avoid
379 # accumulating rounding error
380 work[slice1] = b - matvec(x)
381 resid, info = _stoptest(work[slice1], atol)
382 ijob = 2
384 if info > 0 and iter_ == maxiter and not (resid <= atol):
385 # info isn't set appropriately otherwise
386 info = iter_
388 return postprocess(x), info
391@set_docstring('Use Conjugate Gradient Squared iteration to solve ``Ax = b``.',
392 'The real-valued N-by-N matrix of the linear system.\n'
393 'Alternatively, ``A`` can be a linear operator which can\n'
394 'produce ``Ax`` using, e.g.,\n'
395 '``scipy.sparse.linalg.LinearOperator``.',
396 footer="""\
397 Examples
398 --------
399 >>> import numpy as np
400 >>> from scipy.sparse import csc_matrix
401 >>> from scipy.sparse.linalg import cgs
402 >>> R = np.array([[4, 2, 0, 1],
403 ... [3, 0, 0, 2],
404 ... [0, 1, 1, 1],
405 ... [0, 2, 1, 0]])
406 >>> A = csc_matrix(R)
407 >>> b = np.array([-1, -0.5, -1, 2])
408 >>> x, exit_code = cgs(A, b)
409 >>> print(exit_code) # 0 indicates successful convergence
410 0
411 >>> np.allclose(A.dot(x), b)
412 True
413 """
414 )
415@non_reentrant()
416def cgs(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None):
417 A, M, x, b, postprocess = make_system(A, M, x0, b)
419 n = len(b)
420 if maxiter is None:
421 maxiter = n*10
423 matvec = A.matvec
424 psolve = M.matvec
425 ltr = _type_conv[x.dtype.char]
426 revcom = getattr(_iterative, ltr + 'cgsrevcom')
428 get_residual = lambda: np.linalg.norm(matvec(x) - b)
429 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'cgs')
430 if atol == 'exit':
431 return postprocess(x), 0
433 resid = atol
434 ndx1 = 1
435 ndx2 = -1
436 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
437 work = _aligned_zeros(7*n,dtype=x.dtype)
438 ijob = 1
439 info = 0
440 ftflag = True
441 iter_ = maxiter
442 while True:
443 olditer = iter_
444 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
445 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
446 if callback is not None and iter_ > olditer:
447 callback(x)
448 slice1 = slice(ndx1-1, ndx1-1+n)
449 slice2 = slice(ndx2-1, ndx2-1+n)
450 if (ijob == -1):
451 if callback is not None:
452 callback(x)
453 break
454 elif (ijob == 1):
455 work[slice2] *= sclr2
456 work[slice2] += sclr1*matvec(work[slice1])
457 elif (ijob == 2):
458 work[slice1] = psolve(work[slice2])
459 elif (ijob == 3):
460 work[slice2] *= sclr2
461 work[slice2] += sclr1*matvec(x)
462 elif (ijob == 4):
463 if ftflag:
464 info = -1
465 ftflag = False
466 resid, info = _stoptest(work[slice1], atol)
467 if info == 1 and iter_ > 1:
468 # recompute residual and recheck, to avoid
469 # accumulating rounding error
470 work[slice1] = b - matvec(x)
471 resid, info = _stoptest(work[slice1], atol)
472 ijob = 2
474 if info == -10:
475 # termination due to breakdown: check for convergence
476 resid, ok = _stoptest(b - matvec(x), atol)
477 if ok:
478 info = 0
480 if info > 0 and iter_ == maxiter and not (resid <= atol):
481 # info isn't set appropriately otherwise
482 info = iter_
484 return postprocess(x), info
487@non_reentrant()
488def gmres(A, b, x0=None, tol=1e-5, restart=None, maxiter=None, M=None, callback=None,
489 restrt=None, atol=None, callback_type=None):
490 """
491 Use Generalized Minimal RESidual iteration to solve ``Ax = b``.
493 Parameters
494 ----------
495 A : {sparse matrix, ndarray, LinearOperator}
496 The real or complex N-by-N matrix of the linear system.
497 Alternatively, ``A`` can be a linear operator which can
498 produce ``Ax`` using, e.g.,
499 ``scipy.sparse.linalg.LinearOperator``.
500 b : ndarray
501 Right hand side of the linear system. Has shape (N,) or (N,1).
503 Returns
504 -------
505 x : ndarray
506 The converged solution.
507 info : int
508 Provides convergence information:
509 * 0 : successful exit
510 * >0 : convergence to tolerance not achieved, number of iterations
511 * <0 : illegal input or breakdown
513 Other parameters
514 ----------------
515 x0 : ndarray
516 Starting guess for the solution (a vector of zeros by default).
517 tol, atol : float, optional
518 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``.
519 The default for ``atol`` is ``'legacy'``, which emulates
520 a different legacy behavior.
522 .. warning::
524 The default value for `atol` will be changed in a future release.
525 For future compatibility, specify `atol` explicitly.
526 restart : int, optional
527 Number of iterations between restarts. Larger values increase
528 iteration cost, but may be necessary for convergence.
529 Default is 20.
530 maxiter : int, optional
531 Maximum number of iterations (restart cycles). Iteration will stop
532 after maxiter steps even if the specified tolerance has not been
533 achieved.
534 M : {sparse matrix, ndarray, LinearOperator}
535 Inverse of the preconditioner of A. M should approximate the
536 inverse of A and be easy to solve for (see Notes). Effective
537 preconditioning dramatically improves the rate of convergence,
538 which implies that fewer iterations are needed to reach a given
539 error tolerance. By default, no preconditioner is used.
540 callback : function
541 User-supplied function to call after each iteration. It is called
542 as `callback(args)`, where `args` are selected by `callback_type`.
543 callback_type : {'x', 'pr_norm', 'legacy'}, optional
544 Callback function argument requested:
545 - ``x``: current iterate (ndarray), called on every restart
546 - ``pr_norm``: relative (preconditioned) residual norm (float),
547 called on every inner iteration
548 - ``legacy`` (default): same as ``pr_norm``, but also changes the
549 meaning of 'maxiter' to count inner iterations instead of restart
550 cycles.
551 restrt : int, optional, deprecated
553 .. deprecated:: 0.11.0
554 `gmres` keyword argument `restrt` is deprecated infavour of
555 `restart` and will be removed in SciPy 1.12.0.
557 See Also
558 --------
559 LinearOperator
561 Notes
562 -----
563 A preconditioner, P, is chosen such that P is close to A but easy to solve
564 for. The preconditioner parameter required by this routine is
565 ``M = P^-1``. The inverse should preferably not be calculated
566 explicitly. Rather, use the following template to produce M::
568 # Construct a linear operator that computes P^-1 @ x.
569 import scipy.sparse.linalg as spla
570 M_x = lambda x: spla.spsolve(P, x)
571 M = spla.LinearOperator((n, n), M_x)
573 Examples
574 --------
575 >>> import numpy as np
576 >>> from scipy.sparse import csc_matrix
577 >>> from scipy.sparse.linalg import gmres
578 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
579 >>> b = np.array([2, 4, -1], dtype=float)
580 >>> x, exitCode = gmres(A, b)
581 >>> print(exitCode) # 0 indicates successful convergence
582 0
583 >>> np.allclose(A.dot(x), b)
584 True
585 """
587 # Change 'restrt' keyword to 'restart'
588 if restrt is None:
589 restrt = restart
590 elif restart is not None:
591 raise ValueError("Cannot specify both restart and restrt keywords. "
592 "Preferably use 'restart' only.")
593 else:
594 msg = ("'gmres' keyword argument 'restrt' is deprecated infavour of "
595 "'restart' and will be removed in SciPy 1.12.0.")
596 warnings.warn(msg, DeprecationWarning, stacklevel=2)
598 if callback is not None and callback_type is None:
599 # Warn about 'callback_type' semantic changes.
600 # Probably should be removed only in far future, Scipy 2.0 or so.
601 warnings.warn("scipy.sparse.linalg.gmres called without specifying `callback_type`. "
602 "The default value will be changed in a future release. "
603 "For compatibility, specify a value for `callback_type` explicitly, e.g., "
604 "``{name}(..., callback_type='pr_norm')``, or to retain the old behavior "
605 "``{name}(..., callback_type='legacy')``",
606 category=DeprecationWarning, stacklevel=3)
608 if callback_type is None:
609 callback_type = 'legacy'
611 if callback_type not in ('x', 'pr_norm', 'legacy'):
612 raise ValueError("Unknown callback_type: {!r}".format(callback_type))
614 if callback is None:
615 callback_type = 'none'
617 A, M, x, b,postprocess = make_system(A, M, x0, b)
619 n = len(b)
620 if maxiter is None:
621 maxiter = n*10
623 if restrt is None:
624 restrt = 20
625 restrt = min(restrt, n)
627 matvec = A.matvec
628 psolve = M.matvec
629 ltr = _type_conv[x.dtype.char]
630 revcom = getattr(_iterative, ltr + 'gmresrevcom')
632 bnrm2 = np.linalg.norm(b)
633 Mb_nrm2 = np.linalg.norm(psolve(b))
634 get_residual = lambda: np.linalg.norm(matvec(x) - b)
635 atol = _get_atol(tol, atol, bnrm2, get_residual, 'gmres')
636 if atol == 'exit':
637 return postprocess(x), 0
639 if bnrm2 == 0:
640 return postprocess(b), 0
642 # Tolerance passed to GMRESREVCOM applies to the inner iteration
643 # and deals with the left-preconditioned residual.
644 ptol_max_factor = 1.0
645 ptol = Mb_nrm2 * min(ptol_max_factor, atol / bnrm2)
646 resid = np.nan
647 presid = np.nan
648 ndx1 = 1
649 ndx2 = -1
650 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
651 work = _aligned_zeros((6+restrt)*n,dtype=x.dtype)
652 work2 = _aligned_zeros((restrt+1)*(2*restrt+2),dtype=x.dtype)
653 ijob = 1
654 info = 0
655 ftflag = True
656 iter_ = maxiter
657 old_ijob = ijob
658 first_pass = True
659 resid_ready = False
660 iter_num = 1
661 while True:
662 olditer = iter_
663 x, iter_, presid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
664 revcom(b, x, restrt, work, work2, iter_, presid, info, ndx1, ndx2, ijob, ptol)
665 if callback_type == 'x' and iter_ != olditer:
666 callback(x)
667 slice1 = slice(ndx1-1, ndx1-1+n)
668 slice2 = slice(ndx2-1, ndx2-1+n)
669 if (ijob == -1): # gmres success, update last residual
670 if callback_type in ('pr_norm', 'legacy'):
671 if resid_ready:
672 callback(presid / bnrm2)
673 elif callback_type == 'x':
674 callback(x)
675 break
676 elif (ijob == 1):
677 work[slice2] *= sclr2
678 work[slice2] += sclr1*matvec(x)
679 elif (ijob == 2):
680 work[slice1] = psolve(work[slice2])
681 if not first_pass and old_ijob == 3:
682 resid_ready = True
684 first_pass = False
685 elif (ijob == 3):
686 work[slice2] *= sclr2
687 work[slice2] += sclr1*matvec(work[slice1])
688 if resid_ready:
689 if callback_type in ('pr_norm', 'legacy'):
690 callback(presid / bnrm2)
691 resid_ready = False
692 iter_num = iter_num+1
694 elif (ijob == 4):
695 if ftflag:
696 info = -1
697 ftflag = False
698 resid, info = _stoptest(work[slice1], atol)
700 # Inner loop tolerance control
701 if info or presid > ptol:
702 ptol_max_factor = min(1.0, 1.5 * ptol_max_factor)
703 else:
704 # Inner loop tolerance OK, but outer loop not.
705 ptol_max_factor = max(1e-16, 0.25 * ptol_max_factor)
707 if resid != 0:
708 ptol = presid * min(ptol_max_factor, atol / resid)
709 else:
710 ptol = presid * ptol_max_factor
712 old_ijob = ijob
713 ijob = 2
715 if callback_type == 'legacy':
716 # Legacy behavior
717 if iter_num > maxiter:
718 info = maxiter
719 break
721 if info >= 0 and not (resid <= atol):
722 # info isn't set appropriately otherwise
723 info = maxiter
725 return postprocess(x), info
728@non_reentrant()
729def qmr(A, b, x0=None, tol=1e-5, maxiter=None, M1=None, M2=None, callback=None,
730 atol=None):
731 """Use Quasi-Minimal Residual iteration to solve ``Ax = b``.
733 Parameters
734 ----------
735 A : {sparse matrix, ndarray, LinearOperator}
736 The real-valued N-by-N matrix of the linear system.
737 Alternatively, ``A`` can be a linear operator which can
738 produce ``Ax`` and ``A^T x`` using, e.g.,
739 ``scipy.sparse.linalg.LinearOperator``.
740 b : ndarray
741 Right hand side of the linear system. Has shape (N,) or (N,1).
743 Returns
744 -------
745 x : ndarray
746 The converged solution.
747 info : integer
748 Provides convergence information:
749 0 : successful exit
750 >0 : convergence to tolerance not achieved, number of iterations
751 <0 : illegal input or breakdown
753 Other Parameters
754 ----------------
755 x0 : ndarray
756 Starting guess for the solution.
757 tol, atol : float, optional
758 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``.
759 The default for ``atol`` is ``'legacy'``, which emulates
760 a different legacy behavior.
762 .. warning::
764 The default value for `atol` will be changed in a future release.
765 For future compatibility, specify `atol` explicitly.
766 maxiter : integer
767 Maximum number of iterations. Iteration will stop after maxiter
768 steps even if the specified tolerance has not been achieved.
769 M1 : {sparse matrix, ndarray, LinearOperator}
770 Left preconditioner for A.
771 M2 : {sparse matrix, ndarray, LinearOperator}
772 Right preconditioner for A. Used together with the left
773 preconditioner M1. The matrix M1@A@M2 should have better
774 conditioned than A alone.
775 callback : function
776 User-supplied function to call after each iteration. It is called
777 as callback(xk), where xk is the current solution vector.
779 See Also
780 --------
781 LinearOperator
783 Examples
784 --------
785 >>> import numpy as np
786 >>> from scipy.sparse import csc_matrix
787 >>> from scipy.sparse.linalg import qmr
788 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
789 >>> b = np.array([2, 4, -1], dtype=float)
790 >>> x, exitCode = qmr(A, b)
791 >>> print(exitCode) # 0 indicates successful convergence
792 0
793 >>> np.allclose(A.dot(x), b)
794 True
795 """
796 A_ = A
797 A, M, x, b, postprocess = make_system(A, None, x0, b)
799 if M1 is None and M2 is None:
800 if hasattr(A_,'psolve'):
801 def left_psolve(b):
802 return A_.psolve(b,'left')
804 def right_psolve(b):
805 return A_.psolve(b,'right')
807 def left_rpsolve(b):
808 return A_.rpsolve(b,'left')
810 def right_rpsolve(b):
811 return A_.rpsolve(b,'right')
812 M1 = LinearOperator(A.shape, matvec=left_psolve, rmatvec=left_rpsolve)
813 M2 = LinearOperator(A.shape, matvec=right_psolve, rmatvec=right_rpsolve)
814 else:
815 def id(b):
816 return b
817 M1 = LinearOperator(A.shape, matvec=id, rmatvec=id)
818 M2 = LinearOperator(A.shape, matvec=id, rmatvec=id)
820 n = len(b)
821 if maxiter is None:
822 maxiter = n*10
824 ltr = _type_conv[x.dtype.char]
825 revcom = getattr(_iterative, ltr + 'qmrrevcom')
827 get_residual = lambda: np.linalg.norm(A.matvec(x) - b)
828 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'qmr')
829 if atol == 'exit':
830 return postprocess(x), 0
832 resid = atol
833 ndx1 = 1
834 ndx2 = -1
835 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1
836 work = _aligned_zeros(11*n,x.dtype)
837 ijob = 1
838 info = 0
839 ftflag = True
840 iter_ = maxiter
841 while True:
842 olditer = iter_
843 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
844 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
845 if callback is not None and iter_ > olditer:
846 callback(x)
847 slice1 = slice(ndx1-1, ndx1-1+n)
848 slice2 = slice(ndx2-1, ndx2-1+n)
849 if (ijob == -1):
850 if callback is not None:
851 callback(x)
852 break
853 elif (ijob == 1):
854 work[slice2] *= sclr2
855 work[slice2] += sclr1*A.matvec(work[slice1])
856 elif (ijob == 2):
857 work[slice2] *= sclr2
858 work[slice2] += sclr1*A.rmatvec(work[slice1])
859 elif (ijob == 3):
860 work[slice1] = M1.matvec(work[slice2])
861 elif (ijob == 4):
862 work[slice1] = M2.matvec(work[slice2])
863 elif (ijob == 5):
864 work[slice1] = M1.rmatvec(work[slice2])
865 elif (ijob == 6):
866 work[slice1] = M2.rmatvec(work[slice2])
867 elif (ijob == 7):
868 work[slice2] *= sclr2
869 work[slice2] += sclr1*A.matvec(x)
870 elif (ijob == 8):
871 if ftflag:
872 info = -1
873 ftflag = False
874 resid, info = _stoptest(work[slice1], atol)
875 ijob = 2
877 if info > 0 and iter_ == maxiter and not (resid <= atol):
878 # info isn't set appropriately otherwise
879 info = iter_
881 return postprocess(x), info