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