Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/

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, qr, solve, svd, qr_insert, lstsq)

8from scipy.sparse.linalg._isolve.utils import make_system

11__all__ = ['gcrotmk']

14def _fgmres(matvec, v0, m, atol, lpsolve=None, rpsolve=None, cs=(), outer_v=(),

15 prepend_outer_v=False):

16 """

17 FGMRES Arnoldi process, with optional projection or augmentation

19 Parameters

20 ----------

21 matvec : callable

22 Operation A*x

23 v0 : ndarray

24 Initial vector, normalized to nrm2(v0) == 1

25 m : int

26 Number of GMRES rounds

27 atol : float

28 Absolute tolerance for early exit

29 lpsolve : callable

30 Left preconditioner L

31 rpsolve : callable

32 Right preconditioner R

33 cs : list of (ndarray, ndarray)

34 Columns of matrices C and U in GCROT

35 outer_v : list of ndarrays

36 Augmentation vectors in LGMRES

37 prepend_outer_v : bool, optional

38 Whether augmentation vectors come before or after

39 Krylov iterates

41 Raises

42 ------

43 LinAlgError

44 If nans encountered

46 Returns

47 -------

48 Q, R : ndarray

49 QR decomposition of the upper Hessenberg H=QR

50 B : ndarray

51 Projections corresponding to matrix C

52 vs : list of ndarray

53 Columns of matrix V

54 zs : list of ndarray

55 Columns of matrix Z

56 y : ndarray

57 Solution to ||H y - e_1||_2 = min!

58 res : float

59 The final (preconditioned) residual norm

61 """

63 if lpsolve is None:

64 lpsolve = lambda x: x

65 if rpsolve is None:

66 rpsolve = lambda x: x

68 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (v0,))

70 vs = [v0]

71 zs = []

72 y = None

73 res = np.nan

75 m = m + len(outer_v)

77 # Orthogonal projection coefficients

78 B = np.zeros((len(cs), m), dtype=v0.dtype)

80 # H is stored in QR factorized form

81 Q = np.ones((1, 1), dtype=v0.dtype)

82 R = np.zeros((1, 0), dtype=v0.dtype)

84 eps = np.finfo(v0.dtype).eps

86 breakdown = False

88 # FGMRES Arnoldi process

89 for j in range(m):

90 # L A Z = C B + V H

92 if prepend_outer_v and j < len(outer_v):

93 z, w = outer_v[j]

94 elif prepend_outer_v and j == len(outer_v):

95 z = rpsolve(v0)

96 w = None

97 elif not prepend_outer_v and j >= m - len(outer_v):

98 z, w = outer_v[j - (m - len(outer_v))]

99 else:

100 z = rpsolve(vs[-1])

101 w = None

102

103 if w is None:

104 w = lpsolve(matvec(z))

105 else:

106 # w is clobbered below

107 w = w.copy()

108

109 w_norm = nrm2(w)

110

111 # GCROT projection: L A -> (1 - C C^H) L A

112 # i.e. orthogonalize against C

113 for i, c in enumerate(cs):

114 alpha = dot(c, w)

115 B[i,j] = alpha

116 w = axpy(c, w, c.shape[0], -alpha) # w -= alpha*c

117

118 # Orthogonalize against V

119 hcur = np.zeros(j+2, dtype=Q.dtype)

120 for i, v in enumerate(vs):

121 alpha = dot(v, w)

122 hcur[i] = alpha

123 w = axpy(v, w, v.shape[0], -alpha) # w -= alpha*v

124 hcur[i+1] = nrm2(w)

125

126 with np.errstate(over='ignore', divide='ignore'):

127 # Careful with denormals

128 alpha = 1/hcur[-1]

129

130 if np.isfinite(alpha):

131 w = scal(alpha, w)

132

133 if not (hcur[-1] > eps * w_norm):

134 # w essentially in the span of previous vectors,

135 # or we have nans. Bail out after updating the QR

136 # solution.

137 breakdown = True

138

139 vs.append(w)

140 zs.append(z)

141

142 # Arnoldi LSQ problem

143

144 # Add new column to H=Q@R, padding other columns with zeros

145 Q2 = np.zeros((j+2, j+2), dtype=Q.dtype, order='F')

146 Q2[:j+1,:j+1] = Q

147 Q2[j+1,j+1] = 1

148

149 R2 = np.zeros((j+2, j), dtype=R.dtype, order='F')

150 R2[:j+1,:] = R

151

152 Q, R = qr_insert(Q2, R2, hcur, j, which='col',

153 overwrite_qru=True, check_finite=False)

154

155 # Transformed least squares problem

156 # || Q R y - inner_res_0 * e_1 ||_2 = min!

157 # Since R = [R'; 0], solution is y = inner_res_0 (R')^{-1} (Q^H)[:j,0]

158

159 # Residual is immediately known

160 res = abs(Q[0,-1])

161

162 # Check for termination

163 if res < atol or breakdown:

164 break

165

166 if not np.isfinite(R[j,j]):

167 # nans encountered, bail out

168 raise LinAlgError()

169

170 # -- Get the LSQ problem solution

171

172 # The problem is triangular, but the condition number may be

173 # bad (or in case of breakdown the last diagonal entry may be

174 # zero), so use lstsq instead of trtrs.

175 y, _, _, _, = lstsq(R[:j+1,:j+1], Q[0,:j+1].conj())

176

177 B = B[:,:j+1]

178

179 return Q, R, B, vs, zs, y, res

180

181

182def gcrotmk(A, b, x0=None, tol=1e-5, maxiter=1000, M=None, callback=None,

183 m=20, k=None, CU=None, discard_C=False, truncate='oldest',

184 atol=None):

185 """

186 Solve a matrix equation using flexible GCROT(m,k) algorithm.

187

188 Parameters

189 ----------

190 A : {sparse matrix, ndarray, LinearOperator}

191 The real or complex N-by-N matrix of the linear system.

192 Alternatively, ``A`` can be a linear operator which can

193 produce ``Ax`` using, e.g.,

194 ``scipy.sparse.linalg.LinearOperator``.

195 b : ndarray

196 Right hand side of the linear system. Has shape (N,) or (N,1).

197 x0 : ndarray

198 Starting guess for the solution.

199 tol, atol : float, optional

200 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``.

201 The default for ``atol`` is `tol`.

202

203 .. warning::

204

205 The default value for `atol` will be changed in a future release.

206 For future compatibility, specify `atol` explicitly.

207 maxiter : int, optional

208 Maximum number of iterations. Iteration will stop after maxiter

209 steps even if the specified tolerance has not been achieved.

210 M : {sparse matrix, ndarray, LinearOperator}, optional

211 Preconditioner for A. The preconditioner should approximate the

212 inverse of A. gcrotmk is a 'flexible' algorithm and the preconditioner

213 can vary from iteration to iteration. Effective preconditioning

214 dramatically improves the rate of convergence, which implies that

215 fewer iterations are needed to reach a given error tolerance.

216 callback : function, optional

217 User-supplied function to call after each iteration. It is called

218 as callback(xk), where xk is the current solution vector.

219 m : int, optional

220 Number of inner FGMRES iterations per each outer iteration.

221 Default: 20

222 k : int, optional

223 Number of vectors to carry between inner FGMRES iterations.

224 According to [2]_, good values are around m.

225 Default: m

226 CU : list of tuples, optional

227 List of tuples ``(c, u)`` which contain the columns of the matrices

228 C and U in the GCROT(m,k) algorithm. For details, see [2]_.

229 The list given and vectors contained in it are modified in-place.

230 If not given, start from empty matrices. The ``c`` elements in the

231 tuples can be ``None``, in which case the vectors are recomputed

232 via ``c = A u`` on start and orthogonalized as described in [3]_.

233 discard_C : bool, optional

234 Discard the C-vectors at the end. Useful if recycling Krylov subspaces

235 for different linear systems.

236 truncate : {'oldest', 'smallest'}, optional

237 Truncation scheme to use. Drop: oldest vectors, or vectors with

238 smallest singular values using the scheme discussed in [1,2].

239 See [2]_ for detailed comparison.

240 Default: 'oldest'

241

242 Returns

243 -------

244 x : ndarray

245 The solution found.

246 info : int

247 Provides convergence information:

248

249 * 0 : successful exit

250 * >0 : convergence to tolerance not achieved, number of iterations

251

252 Examples

253 --------

254 >>> import numpy as np

255 >>> from scipy.sparse import csc_matrix

256 >>> from scipy.sparse.linalg import gcrotmk

257 >>> R = np.random.randn(5, 5)

258 >>> A = csc_matrix(R)

259 >>> b = np.random.randn(5)

260 >>> x, exit_code = gcrotmk(A, b)

261 >>> print(exit_code)

262 0

263 >>> np.allclose(A.dot(x), b)

264 True

265

266 References

267 ----------

268 .. [1] E. de Sturler, ''Truncation strategies for optimal Krylov subspace

269 methods'', SIAM J. Numer. Anal. 36, 864 (1999).

270 .. [2] J.E. Hicken and D.W. Zingg, ''A simplified and flexible variant

271 of GCROT for solving nonsymmetric linear systems'',

272 SIAM J. Sci. Comput. 32, 172 (2010).

273 .. [3] M.L. Parks, E. de Sturler, G. Mackey, D.D. Johnson, S. Maiti,

274 ''Recycling Krylov subspaces for sequences of linear systems'',

275 SIAM J. Sci. Comput. 28, 1651 (2006).

276

277 """

278 A,M,x,b,postprocess = make_system(A,M,x0,b)

279

280 if not np.isfinite(b).all():

281 raise ValueError("RHS must contain only finite numbers")

282

283 if truncate not in ('oldest', 'smallest'):

284 raise ValueError("Invalid value for 'truncate': %r" % (truncate,))

285

286 if atol is None:

287 warnings.warn("scipy.sparse.linalg.gcrotmk called without specifying `atol`. "

288 "The default value will change in the future. To preserve "

289 "current behavior, set ``atol=tol``.",

290 category=DeprecationWarning, stacklevel=2)

291 atol = tol

292

293 matvec = A.matvec

294 psolve = M.matvec

295

296 if CU is None:

297 CU = []

298

299 if k is None:

300 k = m

301

302 axpy, dot, scal = None, None, None

303

304 if x0 is None:

305 r = b.copy()

306 else:

307 r = b - matvec(x)

308

309 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (x, r))

310

311 b_norm = nrm2(b)

312 if b_norm == 0:

313 x = b

314 return (postprocess(x), 0)

315

316 if discard_C:

317 CU[:] = [(None, u) for c, u in CU]

318

319 # Reorthogonalize old vectors

320 if CU:

321 # Sort already existing vectors to the front

322 CU.sort(key=lambda cu: cu[0] is not None)

323

324 # Fill-in missing ones

325 C = np.empty((A.shape[0], len(CU)), dtype=r.dtype, order='F')

326 us = []

327 j = 0

328 while CU:

329 # More memory-efficient: throw away old vectors as we go

330 c, u = CU.pop(0)

331 if c is None:

332 c = matvec(u)

333 C[:,j] = c

334 j += 1

335 us.append(u)

336

337 # Orthogonalize

338 Q, R, P = qr(C, overwrite_a=True, mode='economic', pivoting=True)

339 del C

340

341 # C := Q

342 cs = list(Q.T)

343

344 # U := U P R^-1, back-substitution

345 new_us = []

346 for j in range(len(cs)):

347 u = us[P[j]]

348 for i in range(j):

349 u = axpy(us[P[i]], u, u.shape[0], -R[i,j])

350 if abs(R[j,j]) < 1e-12 * abs(R[0,0]):

351 # discard rest of the vectors

352 break

353 u = scal(1.0/R[j,j], u)

354 new_us.append(u)

355

356 # Form the new CU lists

357 CU[:] = list(zip(cs, new_us))[::-1]

358

359 if CU:

360 axpy, dot = get_blas_funcs(['axpy', 'dot'], (r,))

361

362 # Solve first the projection operation with respect to the CU

363 # vectors. This corresponds to modifying the initial guess to

364 # be

365 #

366 # x' = x + U y

367 # y = argmin_y || b - A (x + U y) ||^2

368 #

369 # The solution is y = C^H (b - A x)

370 for c, u in CU:

371 yc = dot(c, r)

372 x = axpy(u, x, x.shape[0], yc)

373 r = axpy(c, r, r.shape[0], -yc)

374

375 # GCROT main iteration

376 for j_outer in range(maxiter):

377 # -- callback

378 if callback is not None:

379 callback(x)

380

381 beta = nrm2(r)

382

383 # -- check stopping condition

384 beta_tol = max(atol, tol * b_norm)

385

386 if beta <= beta_tol and (j_outer > 0 or CU):

387 # recompute residual to avoid rounding error

388 r = b - matvec(x)

389 beta = nrm2(r)

390

391 if beta <= beta_tol:

392 j_outer = -1

393 break

394

395 ml = m + max(k - len(CU), 0)

396

397 cs = [c for c, u in CU]

398

399 try:

400 Q, R, B, vs, zs, y, pres = _fgmres(matvec,

401 r/beta,

402 ml,

403 rpsolve=psolve,

404 atol=max(atol, tol*b_norm)/beta,

405 cs=cs)

406 y *= beta

407 except LinAlgError:

408 # Floating point over/underflow, non-finite result from

409 # matmul etc. -- report failure.

410 break

411

412 #

413 # At this point,

414 #

415 # [A U, A Z] = [C, V] G; G = [ I B ]

416 # [ 0 H ]

417 #

418 # where [C, V] has orthonormal columns, and r = beta v_0. Moreover,

419 #

420 # || b - A (x + Z y + U q) ||_2 = || r - C B y - V H y - C q ||_2 = min!

421 #

422 # from which y = argmin_y || beta e_1 - H y ||_2, and q = -B y

423 #

424

425 #

426 # GCROT(m,k) update

427 #

428

429 # Define new outer vectors

430

431 # ux := (Z - U B) y

432 ux = zs[0]*y[0]

433 for z, yc in zip(zs[1:], y[1:]):

434 ux = axpy(z, ux, ux.shape[0], yc) # ux += z*yc

435 by = B.dot(y)

436 for cu, byc in zip(CU, by):

437 c, u = cu

438 ux = axpy(u, ux, ux.shape[0], -byc) # ux -= u*byc

439

440 # cx := V H y

441 hy = Q.dot(R.dot(y))

442 cx = vs[0] * hy[0]

443 for v, hyc in zip(vs[1:], hy[1:]):

444 cx = axpy(v, cx, cx.shape[0], hyc) # cx += v*hyc

445

446 # Normalize cx, maintaining cx = A ux

447 # This new cx is orthogonal to the previous C, by construction

448 try:

449 alpha = 1/nrm2(cx)

450 if not np.isfinite(alpha):

451 raise FloatingPointError()

452 except (FloatingPointError, ZeroDivisionError):

453 # Cannot update, so skip it

454 continue

455

456 cx = scal(alpha, cx)

457 ux = scal(alpha, ux)

458

459 # Update residual and solution

460 gamma = dot(cx, r)

461 r = axpy(cx, r, r.shape[0], -gamma) # r -= gamma*cx

462 x = axpy(ux, x, x.shape[0], gamma) # x += gamma*ux

463

464 # Truncate CU

465 if truncate == 'oldest':

466 while len(CU) >= k and CU:

467 del CU[0]

468 elif truncate == 'smallest':

469 if len(CU) >= k and CU:

470 # cf. [1,2]

471 D = solve(R[:-1,:].T, B.T).T

472 W, sigma, V = svd(D)

473

474 # C := C W[:,:k-1], U := U W[:,:k-1]

475 new_CU = []

476 for j, w in enumerate(W[:,:k-1].T):

477 c, u = CU[0]

478 c = c * w[0]

479 u = u * w[0]

480 for cup, wp in zip(CU[1:], w[1:]):

481 cp, up = cup

482 c = axpy(cp, c, c.shape[0], wp)

483 u = axpy(up, u, u.shape[0], wp)

484

485 # Reorthogonalize at the same time; not necessary

486 # in exact arithmetic, but floating point error

487 # tends to accumulate here

488 for cp, up in new_CU:

489 alpha = dot(cp, c)

490 c = axpy(cp, c, c.shape[0], -alpha)

491 u = axpy(up, u, u.shape[0], -alpha)

492 alpha = nrm2(c)

493 c = scal(1.0/alpha, c)

494 u = scal(1.0/alpha, u)

495

496 new_CU.append((c, u))

497 CU[:] = new_CU

498

499 # Add new vector to CU

500 CU.append((cx, ux))

501

502 # Include the solution vector to the span

503 CU.append((None, x.copy()))

504 if discard_C:

505 CU[:] = [(None, uz) for cz, uz in CU]

506

507 return postprocess(x), j_outer + 1

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/_gcrotmk.py: 4%

195 statements