Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/sparse/linalg/isolve/lgmres.py : 12%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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, dense matrix, 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 : {array, matrix}
33 Right hand side of the linear system. Has shape (N,) or (N,1).
34 x0 : {array, matrix}
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, dense matrix, 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 : array or matrix
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 >>> from scipy.sparse import csc_matrix
116 >>> from scipy.sparse.linalg import lgmres
117 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
118 >>> b = np.array([2, 4, -1], dtype=float)
119 >>> x, exitCode = lgmres(A, b)
120 >>> print(exitCode) # 0 indicates successful convergence
121 0
122 >>> np.allclose(A.dot(x), b)
123 True
124 """
125 A,M,x,b,postprocess = make_system(A,M,x0,b)
127 if not np.isfinite(b).all():
128 raise ValueError("RHS must contain only finite numbers")
130 if atol is None:
131 warnings.warn("scipy.sparse.linalg.lgmres called without specifying `atol`. "
132 "The default value will change in the future. To preserve "
133 "current behavior, set ``atol=tol``.",
134 category=DeprecationWarning, stacklevel=2)
135 atol = tol
137 matvec = A.matvec
138 psolve = M.matvec
140 if outer_v is None:
141 outer_v = []
143 axpy, dot, scal = None, None, None
144 nrm2 = get_blas_funcs('nrm2', [b])
146 b_norm = nrm2(b)
147 ptol_max_factor = 1.0
149 for k_outer in range(maxiter):
150 r_outer = matvec(x) - b
152 # -- callback
153 if callback is not None:
154 callback(x)
156 # -- determine input type routines
157 if axpy is None:
158 if np.iscomplexobj(r_outer) and not np.iscomplexobj(x):
159 x = x.astype(r_outer.dtype)
160 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'],
161 (x, r_outer))
163 # -- check stopping condition
164 r_norm = nrm2(r_outer)
165 if r_norm <= max(atol, tol * b_norm):
166 break
168 # -- inner LGMRES iteration
169 v0 = -psolve(r_outer)
170 inner_res_0 = nrm2(v0)
172 if inner_res_0 == 0:
173 rnorm = nrm2(r_outer)
174 raise RuntimeError("Preconditioner returned a zero vector; "
175 "|v| ~ %.1g, |M v| = 0" % rnorm)
177 v0 = scal(1.0/inner_res_0, v0)
179 ptol = min(ptol_max_factor, max(atol, tol*b_norm)/r_norm)
181 try:
182 Q, R, B, vs, zs, y, pres = _fgmres(matvec,
183 v0,
184 inner_m,
185 lpsolve=psolve,
186 atol=ptol,
187 outer_v=outer_v,
188 prepend_outer_v=prepend_outer_v)
189 y *= inner_res_0
190 if not np.isfinite(y).all():
191 # Overflow etc. in computation. There's no way to
192 # recover from this, so we have to bail out.
193 raise LinAlgError()
194 except LinAlgError:
195 # Floating point over/underflow, non-finite result from
196 # matmul etc. -- report failure.
197 return postprocess(x), k_outer + 1
199 # Inner loop tolerance control
200 if pres > ptol:
201 ptol_max_factor = min(1.0, 1.5 * ptol_max_factor)
202 else:
203 ptol_max_factor = max(1e-16, 0.25 * ptol_max_factor)
205 # -- GMRES terminated: eval solution
206 dx = zs[0]*y[0]
207 for w, yc in zip(zs[1:], y[1:]):
208 dx = axpy(w, dx, dx.shape[0], yc) # dx += w*yc
210 # -- Store LGMRES augmentation vectors
211 nx = nrm2(dx)
212 if nx > 0:
213 if store_outer_Av:
214 q = Q.dot(R.dot(y))
215 ax = vs[0]*q[0]
216 for v, qc in zip(vs[1:], q[1:]):
217 ax = axpy(v, ax, ax.shape[0], qc)
218 outer_v.append((dx/nx, ax/nx))
219 else:
220 outer_v.append((dx/nx, None))
222 # -- Retain only a finite number of augmentation vectors
223 while len(outer_v) > outer_k:
224 del outer_v[0]
226 # -- Apply step
227 x += dx
228 else:
229 # didn't converge ...
230 return postprocess(x), maxiter
232 return postprocess(x), 0