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

2# Distributed under the same license as SciPy.

4import numpy as np

5from numpy.linalg import LinAlgError

6from scipy.linalg import (get_blas_funcs, qr, solve, svd, qr_insert, lstsq)

7from .iterative import _get_atol_rtol

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

9from scipy._lib.deprecation import _NoValue, _deprecate_positional_args

12__all__ = ['gcrotmk']

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

16 prepend_outer_v=False):

17 """

18 FGMRES Arnoldi process, with optional projection or augmentation

20 Parameters

21 ----------

22 matvec : callable

23 Operation A*x

24 v0 : ndarray

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

26 m : int

27 Number of GMRES rounds

28 atol : float

29 Absolute tolerance for early exit

30 lpsolve : callable

31 Left preconditioner L

32 rpsolve : callable

33 Right preconditioner R

34 cs : list of (ndarray, ndarray)

35 Columns of matrices C and U in GCROT

36 outer_v : list of ndarrays

37 Augmentation vectors in LGMRES

38 prepend_outer_v : bool, optional

39 Whether augmentation vectors come before or after

40 Krylov iterates

42 Raises

43 ------

44 LinAlgError

45 If nans encountered

47 Returns

48 -------

49 Q, R : ndarray

50 QR decomposition of the upper Hessenberg H=QR

51 B : ndarray

52 Projections corresponding to matrix C

53 vs : list of ndarray

54 Columns of matrix V

55 zs : list of ndarray

56 Columns of matrix Z

57 y : ndarray

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

59 res : float

60 The final (preconditioned) residual norm

62 """

64 if lpsolve is None:

65 def lpsolve(x):

66 return x

67 if rpsolve is None:

68 def rpsolve(x):

69 return x

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

73 vs = [v0]

74 zs = []

75 y = None

76 res = np.nan

78 m = m + len(outer_v)

80 # Orthogonal projection coefficients

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

83 # H is stored in QR factorized form

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

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

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

89 breakdown = False

91 # FGMRES Arnoldi process

92 for j in range(m):

93 # L A Z = C B + V H

95 if prepend_outer_v and j < len(outer_v):

96 z, w = outer_v[j]

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

98 z = rpsolve(v0)

99 w = None

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

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

102 else:

103 z = rpsolve(vs[-1])

104 w = None

105

106 if w is None:

107 w = lpsolve(matvec(z))

108 else:

109 # w is clobbered below

110 w = w.copy()

111

112 w_norm = nrm2(w)

113

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

115 # i.e. orthogonalize against C

116 for i, c in enumerate(cs):

117 alpha = dot(c, w)

118 B[i,j] = alpha

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

120

121 # Orthogonalize against V

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

123 for i, v in enumerate(vs):

124 alpha = dot(v, w)

125 hcur[i] = alpha

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

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

128

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

130 # Careful with denormals

131 alpha = 1/hcur[-1]

132

133 if np.isfinite(alpha):

134 w = scal(alpha, w)

135

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

137 # w essentially in the span of previous vectors,

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

139 # solution.

140 breakdown = True

141

142 vs.append(w)

143 zs.append(z)

144

145 # Arnoldi LSQ problem

146

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

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

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

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

151

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

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

154

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

156 overwrite_qru=True, check_finite=False)

157

158 # Transformed least squares problem

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

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

161

162 # Residual is immediately known

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

164

165 # Check for termination

166 if res < atol or breakdown:

167 break

168

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

170 # nans encountered, bail out

171 raise LinAlgError()

172

173 # -- Get the LSQ problem solution

174

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

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

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

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

179

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

181

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

183

184

185@_deprecate_positional_args(version="1.14.0")

186def gcrotmk(A, b, x0=None, *, tol=_NoValue, maxiter=1000, M=None, callback=None,

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

188 atol=None, rtol=1e-5):

189 """

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

191

192 Parameters

193 ----------

194 A : {sparse matrix, ndarray, LinearOperator}

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

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

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

198 ``scipy.sparse.linalg.LinearOperator``.

199 b : ndarray

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

201 x0 : ndarray

202 Starting guess for the solution.

203 rtol, atol : float, optional

204 Parameters for the convergence test. For convergence,

205 ``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.

206 The default is ``rtol=1e-5``, the default for ``atol`` is ``rtol``.

207

208 .. warning::

209

210 The default value for ``atol`` will be changed to ``0.0`` in

211 SciPy 1.14.0.

212 maxiter : int, optional

213 Maximum number of iterations. Iteration will stop after maxiter

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

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

216 Preconditioner for A. The preconditioner should approximate the

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

218 can vary from iteration to iteration. Effective preconditioning

219 dramatically improves the rate of convergence, which implies that

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

221 callback : function, optional

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

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

224 m : int, optional

225 Number of inner FGMRES iterations per each outer iteration.

226 Default: 20

227 k : int, optional

228 Number of vectors to carry between inner FGMRES iterations.

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

230 Default: m

231 CU : list of tuples, optional

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

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

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

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

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

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

238 discard_C : bool, optional

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

240 for different linear systems.

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

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

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

244 See [2]_ for detailed comparison.

245 Default: 'oldest'

246 tol : float, optional, deprecated

247

248 .. deprecated:: 1.12.0

249 `gcrotmk` keyword argument ``tol`` is deprecated in favor of

250 ``rtol`` and will be removed in SciPy 1.14.0.

251

252 Returns

253 -------

254 x : ndarray

255 The solution found.

256 info : int

257 Provides convergence information:

258

259 * 0 : successful exit

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

261

262 Examples

263 --------

264 >>> import numpy as np

265 >>> from scipy.sparse import csc_matrix

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

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

268 >>> A = csc_matrix(R)

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

270 >>> x, exit_code = gcrotmk(A, b, atol=1e-5)

271 >>> print(exit_code)

272 0

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

274 True

275

276 References

277 ----------

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

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

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

281 of GCROT for solving nonsymmetric linear systems'',

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

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

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

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

286

287 """

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

289

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

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

292

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

294 raise ValueError(f"Invalid value for 'truncate': {truncate!r}")

295

296 matvec = A.matvec

297 psolve = M.matvec

298

299 if CU is None:

300 CU = []

301

302 if k is None:

303 k = m

304

305 axpy, dot, scal = None, None, None

306

307 if x0 is None:

308 r = b.copy()

309 else:

310 r = b - matvec(x)

311

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

313

314 b_norm = nrm2(b)

315

316 # we call this to get the right atol/rtol and raise warnings as necessary

317 atol, rtol = _get_atol_rtol('gcrotmk', b_norm, tol, atol, rtol)

318

319 if b_norm == 0:

320 x = b

321 return (postprocess(x), 0)

322

323 if discard_C:

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

325

326 # Reorthogonalize old vectors

327 if CU:

328 # Sort already existing vectors to the front

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

330

331 # Fill-in missing ones

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

333 us = []

334 j = 0

335 while CU:

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

337 c, u = CU.pop(0)

338 if c is None:

339 c = matvec(u)

340 C[:,j] = c

341 j += 1

342 us.append(u)

343

344 # Orthogonalize

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

346 del C

347

348 # C := Q

349 cs = list(Q.T)

350

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

352 new_us = []

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

354 u = us[P[j]]

355 for i in range(j):

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

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

358 # discard rest of the vectors

359 break

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

361 new_us.append(u)

362

363 # Form the new CU lists

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

365

366 if CU:

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

368

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

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

371 # be

372 #

373 # x' = x + U y

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

375 #

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

377 for c, u in CU:

378 yc = dot(c, r)

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

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

381

382 # GCROT main iteration

383 for j_outer in range(maxiter):

384 # -- callback

385 if callback is not None:

386 callback(x)

387

388 beta = nrm2(r)

389

390 # -- check stopping condition

391 beta_tol = max(atol, rtol * b_norm)

392

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

394 # recompute residual to avoid rounding error

395 r = b - matvec(x)

396 beta = nrm2(r)

397

398 if beta <= beta_tol:

399 j_outer = -1

400 break

401

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

403

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

405

406 try:

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

408 r/beta,

409 ml,

410 rpsolve=psolve,

411 atol=max(atol, rtol*b_norm)/beta,

412 cs=cs)

413 y *= beta

414 except LinAlgError:

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

416 # matmul etc. -- report failure.

417 break

418

419 #

420 # At this point,

421 #

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

423 # [ 0 H ]

424 #

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

426 #

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

428 #

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

430 #

431

432 #

433 # GCROT(m,k) update

434 #

435

436 # Define new outer vectors

437

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

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

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

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

442 by = B.dot(y)

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

444 c, u = cu

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

446

447 # cx := V H y

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

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

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

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

452

453 # Normalize cx, maintaining cx = A ux

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

455 try:

456 alpha = 1/nrm2(cx)

457 if not np.isfinite(alpha):

458 raise FloatingPointError()

459 except (FloatingPointError, ZeroDivisionError):

460 # Cannot update, so skip it

461 continue

462

463 cx = scal(alpha, cx)

464 ux = scal(alpha, ux)

465

466 # Update residual and solution

467 gamma = dot(cx, r)

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

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

470

471 # Truncate CU

472 if truncate == 'oldest':

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

474 del CU[0]

475 elif truncate == 'smallest':

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

477 # cf. [1,2]

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

479 W, sigma, V = svd(D)

480

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

482 new_CU = []

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

484 c, u = CU[0]

485 c = c * w[0]

486 u = u * w[0]

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

488 cp, up = cup

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

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

491

492 # Reorthogonalize at the same time; not necessary

493 # in exact arithmetic, but floating point error

494 # tends to accumulate here

495 for cp, up in new_CU:

496 alpha = dot(cp, c)

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

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

499 alpha = nrm2(c)

500 c = scal(1.0/alpha, c)

501 u = scal(1.0/alpha, u)

502

503 new_CU.append((c, u))

504 CU[:] = new_CU

505

506 # Add new vector to CU

507 CU.append((cx, ux))

508

509 # Include the solution vector to the span

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

511 if discard_C:

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

513

514 return postprocess(x), j_outer + 1

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_isolve/_gcrotmk.py: 5%

197 statements