Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_isolve/lgmres.py: 14%
70 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-03-22 06:44 +0000
1# Copyright (C) 2009, Pauli Virtanen <pav@iki.fi>
2# Distributed under the same license as SciPy.
4import numpy as np
5from numpy.linalg import LinAlgError
6from scipy.linalg import get_blas_funcs
7from .iterative import _get_atol_rtol
8from .utils import make_system
9from scipy._lib.deprecation import _NoValue, _deprecate_positional_args
11from ._gcrotmk import _fgmres
13__all__ = ['lgmres']
16@_deprecate_positional_args(version="1.14.0")
17def lgmres(A, b, x0=None, *, tol=_NoValue, maxiter=1000, M=None, callback=None,
18 inner_m=30, outer_k=3, outer_v=None, store_outer_Av=True,
19 prepend_outer_v=False, atol=None, rtol=1e-5):
20 """
21 Solve a matrix equation using the LGMRES algorithm.
23 The LGMRES algorithm [1]_ [2]_ is designed to avoid some problems
24 in the convergence in restarted GMRES, and often converges in fewer
25 iterations.
27 Parameters
28 ----------
29 A : {sparse matrix, ndarray, LinearOperator}
30 The real or complex N-by-N matrix of the linear system.
31 Alternatively, ``A`` can be a linear operator which can
32 produce ``Ax`` using, e.g.,
33 ``scipy.sparse.linalg.LinearOperator``.
34 b : ndarray
35 Right hand side of the linear system. Has shape (N,) or (N,1).
36 x0 : ndarray
37 Starting guess for the solution.
38 rtol, atol : float, optional
39 Parameters for the convergence test. For convergence,
40 ``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
41 The default is ``rtol=1e-5``, the default for ``atol`` is ``rtol``.
43 .. warning::
45 The default value for ``atol`` will be changed to ``0.0`` in
46 SciPy 1.14.0.
47 maxiter : int, optional
48 Maximum number of iterations. Iteration will stop after maxiter
49 steps even if the specified tolerance has not been achieved.
50 M : {sparse matrix, ndarray, LinearOperator}, optional
51 Preconditioner for A. The preconditioner should approximate the
52 inverse of A. Effective preconditioning dramatically improves the
53 rate of convergence, which implies that fewer iterations are needed
54 to reach a given error tolerance.
55 callback : function, optional
56 User-supplied function to call after each iteration. It is called
57 as callback(xk), where xk is the current solution vector.
58 inner_m : int, optional
59 Number of inner GMRES iterations per each outer iteration.
60 outer_k : int, optional
61 Number of vectors to carry between inner GMRES iterations.
62 According to [1]_, good values are in the range of 1...3.
63 However, note that if you want to use the additional vectors to
64 accelerate solving multiple similar problems, larger values may
65 be beneficial.
66 outer_v : list of tuples, optional
67 List containing tuples ``(v, Av)`` of vectors and corresponding
68 matrix-vector products, used to augment the Krylov subspace, and
69 carried between inner GMRES iterations. The element ``Av`` can
70 be `None` if the matrix-vector product should be re-evaluated.
71 This parameter is modified in-place by `lgmres`, and can be used
72 to pass "guess" vectors in and out of the algorithm when solving
73 similar problems.
74 store_outer_Av : bool, optional
75 Whether LGMRES should store also A@v in addition to vectors `v`
76 in the `outer_v` list. Default is True.
77 prepend_outer_v : bool, optional
78 Whether to put outer_v augmentation vectors before Krylov iterates.
79 In standard LGMRES, prepend_outer_v=False.
80 tol : float, optional, deprecated
82 .. deprecated:: 1.12.0
83 `lgmres` keyword argument ``tol`` is deprecated in favor of ``rtol``
84 and will be removed in SciPy 1.14.0.
86 Returns
87 -------
88 x : ndarray
89 The converged solution.
90 info : int
91 Provides convergence information:
93 - 0 : successful exit
94 - >0 : convergence to tolerance not achieved, number of iterations
95 - <0 : illegal input or breakdown
97 Notes
98 -----
99 The LGMRES algorithm [1]_ [2]_ is designed to avoid the
100 slowing of convergence in restarted GMRES, due to alternating
101 residual vectors. Typically, it often outperforms GMRES(m) of
102 comparable memory requirements by some measure, or at least is not
103 much worse.
105 Another advantage in this algorithm is that you can supply it with
106 'guess' vectors in the `outer_v` argument that augment the Krylov
107 subspace. If the solution lies close to the span of these vectors,
108 the algorithm converges faster. This can be useful if several very
109 similar matrices need to be inverted one after another, such as in
110 Newton-Krylov iteration where the Jacobian matrix often changes
111 little in the nonlinear steps.
113 References
114 ----------
115 .. [1] A.H. Baker and E.R. Jessup and T. Manteuffel, "A Technique for
116 Accelerating the Convergence of Restarted GMRES", SIAM J. Matrix
117 Anal. Appl. 26, 962 (2005).
118 .. [2] A.H. Baker, "On Improving the Performance of the Linear Solver
119 restarted GMRES", PhD thesis, University of Colorado (2003).
121 Examples
122 --------
123 >>> import numpy as np
124 >>> from scipy.sparse import csc_matrix
125 >>> from scipy.sparse.linalg import lgmres
126 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
127 >>> b = np.array([2, 4, -1], dtype=float)
128 >>> x, exitCode = lgmres(A, b, atol=1e-5)
129 >>> print(exitCode) # 0 indicates successful convergence
130 0
131 >>> np.allclose(A.dot(x), b)
132 True
133 """
134 A,M,x,b,postprocess = make_system(A,M,x0,b)
136 if not np.isfinite(b).all():
137 raise ValueError("RHS must contain only finite numbers")
139 matvec = A.matvec
140 psolve = M.matvec
142 if outer_v is None:
143 outer_v = []
145 axpy, dot, scal = None, None, None
146 nrm2 = get_blas_funcs('nrm2', [b])
148 b_norm = nrm2(b)
150 # we call this to get the right atol/rtol and raise warnings as necessary
151 atol, rtol = _get_atol_rtol('lgmres', b_norm, tol, atol, rtol)
153 if b_norm == 0:
154 x = b
155 return (postprocess(x), 0)
157 ptol_max_factor = 1.0
159 for k_outer in range(maxiter):
160 r_outer = matvec(x) - b
162 # -- callback
163 if callback is not None:
164 callback(x)
166 # -- determine input type routines
167 if axpy is None:
168 if np.iscomplexobj(r_outer) and not np.iscomplexobj(x):
169 x = x.astype(r_outer.dtype)
170 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'],
171 (x, r_outer))
173 # -- check stopping condition
174 r_norm = nrm2(r_outer)
175 if r_norm <= max(atol, rtol * b_norm):
176 break
178 # -- inner LGMRES iteration
179 v0 = -psolve(r_outer)
180 inner_res_0 = nrm2(v0)
182 if inner_res_0 == 0:
183 rnorm = nrm2(r_outer)
184 raise RuntimeError("Preconditioner returned a zero vector; "
185 "|v| ~ %.1g, |M v| = 0" % rnorm)
187 v0 = scal(1.0/inner_res_0, v0)
189 ptol = min(ptol_max_factor, max(atol, rtol*b_norm)/r_norm)
191 try:
192 Q, R, B, vs, zs, y, pres = _fgmres(matvec,
193 v0,
194 inner_m,
195 lpsolve=psolve,
196 atol=ptol,
197 outer_v=outer_v,
198 prepend_outer_v=prepend_outer_v)
199 y *= inner_res_0
200 if not np.isfinite(y).all():
201 # Overflow etc. in computation. There's no way to
202 # recover from this, so we have to bail out.
203 raise LinAlgError()
204 except LinAlgError:
205 # Floating point over/underflow, non-finite result from
206 # matmul etc. -- report failure.
207 return postprocess(x), k_outer + 1
209 # Inner loop tolerance control
210 if pres > ptol:
211 ptol_max_factor = min(1.0, 1.5 * ptol_max_factor)
212 else:
213 ptol_max_factor = max(1e-16, 0.25 * ptol_max_factor)
215 # -- GMRES terminated: eval solution
216 dx = zs[0]*y[0]
217 for w, yc in zip(zs[1:], y[1:]):
218 dx = axpy(w, dx, dx.shape[0], yc) # dx += w*yc
220 # -- Store LGMRES augmentation vectors
221 nx = nrm2(dx)
222 if nx > 0:
223 if store_outer_Av:
224 q = Q.dot(R.dot(y))
225 ax = vs[0]*q[0]
226 for v, qc in zip(vs[1:], q[1:]):
227 ax = axpy(v, ax, ax.shape[0], qc)
228 outer_v.append((dx/nx, ax/nx))
229 else:
230 outer_v.append((dx/nx, None))
232 # -- Retain only a finite number of augmentation vectors
233 while len(outer_v) > outer_k:
234 del outer_v[0]
236 # -- Apply step
237 x += dx
238 else:
239 # didn't converge ...
240 return postprocess(x), maxiter
242 return postprocess(x), 0