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

1import warnings 

2import numpy as np 

3from scipy.sparse.linalg._interface import LinearOperator 

4from .utils import make_system 

5from scipy.linalg import get_lapack_funcs 

6from scipy._lib.deprecation import _NoValue, _deprecate_positional_args 

7 

8__all__ = ['bicg', 'bicgstab', 'cg', 'cgs', 'gmres', 'qmr'] 

9 

10 

11def _get_atol(name, b, tol=_NoValue, atol=0., rtol=1e-5): 

12 """ 

13 A helper function to handle tolerance deprecations and normalization 

14 """ 

15 if tol is not _NoValue: 

16 msg = (f"'scipy.sparse.linalg.{name}' keyword argument 'tol' is " 

17 "deprecated in favor of 'rtol' and will be removed in SciPy " 

18 "v.1.14.0. Until then, if set, it will override 'rtol'.") 

19 warnings.warn(msg, category=DeprecationWarning, stacklevel=4) 

20 rtol = float(tol) if tol is not None else rtol 

21 

22 if atol == 'legacy': 

23 warnings.warn("scipy.sparse.linalg.{name} called with `atol` set to " 

24 "string. This behavior is deprecated and atol parameter" 

25 " only excepts floats. In SciPy 1.14, this will result" 

26 " with an error.", category=DeprecationWarning, 

27 stacklevel=4) 

28 atol = 0 

29 

30 atol = max(float(atol), float(rtol) * float(np.linalg.norm(b))) 

31 

32 return atol 

33 

34 

35@_deprecate_positional_args(version="1.14") 

36def bicg(A, b, x0=None, *, tol=_NoValue, maxiter=None, M=None, callback=None, 

37 atol=0., rtol=1e-5): 

38 """Use BIConjugate Gradient iteration to solve ``Ax = b``. 

39 

40 Parameters 

41 ---------- 

42 A : {sparse matrix, ndarray, LinearOperator} 

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

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

45 produce ``Ax`` and ``A^T x`` using, e.g., 

46 ``scipy.sparse.linalg.LinearOperator``. 

47 b : ndarray 

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

49 x0 : ndarray 

50 Starting guess for the solution. 

51 rtol, atol : float, optional 

52 Parameters for the convergence test. For convergence, 

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

54 The default is ``atol=0.`` and ``rtol=1e-5``. 

55 maxiter : integer 

56 Maximum number of iterations. Iteration will stop after maxiter 

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

58 M : {sparse matrix, ndarray, LinearOperator} 

59 Preconditioner for A. The preconditioner should approximate the 

60 inverse of A. Effective preconditioning dramatically improves the 

61 rate of convergence, which implies that fewer iterations are needed 

62 to reach a given error tolerance. 

63 callback : function 

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

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

66 tol : float, optional, deprecated 

67 

68 .. deprecated 1.12.0 

69 `bicg` keyword argument `tol` is deprecated in favor of `rtol` and 

70 will be removed in SciPy 1.14.0. 

71 

72 Returns 

73 ------- 

74 x : ndarray 

75 The converged solution. 

76 info : integer 

77 Provides convergence information: 

78 0 : successful exit 

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

80 <0 : parameter breakdown 

81 

82 Examples 

83 -------- 

84 >>> import numpy as np 

85 >>> from scipy.sparse import csc_matrix 

86 >>> from scipy.sparse.linalg import bicg 

87 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1.]]) 

88 >>> b = np.array([2., 4., -1.]) 

89 >>> x, exitCode = bicg(A, b, atol=1e-5) 

90 >>> print(exitCode) # 0 indicates successful convergence 

91 0 

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

93 True 

94 

95 """ 

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

97 bnrm2 = np.linalg.norm(b) 

98 

99 if bnrm2 == 0: 

100 return postprocess(b), 0 

101 

102 atol = _get_atol('bicg', b, tol, atol, rtol) 

103 

104 n = len(b) 

105 dotprod = np.vdot if np.iscomplexobj(x) else np.dot 

106 

107 if maxiter is None: 

108 maxiter = n*10 

109 

110 matvec, rmatvec = A.matvec, A.rmatvec 

111 psolve, rpsolve = M.matvec, M.rmatvec 

112 

113 rhotol = np.finfo(x.dtype.char).eps**2 

114 

115 # Dummy values to initialize vars, silence linter warnings 

116 rho_prev, p, ptilde = None, None, None 

117 

118 r = b - matvec(x) if x.any() else b.copy() 

119 rtilde = r.copy() 

120 

121 for iteration in range(maxiter): 

122 if np.linalg.norm(r) < atol: # Are we done? 

123 return postprocess(x), 0 

124 

125 z = psolve(r) 

126 ztilde = rpsolve(rtilde) 

127 # order matters in this dot product 

128 rho_cur = dotprod(rtilde, z) 

129 

130 if np.abs(rho_cur) < rhotol: # Breakdown case 

131 return postprocess, -10 

132 

133 if iteration > 0: 

134 beta = rho_cur / rho_prev 

135 p *= beta 

136 p += z 

137 ptilde *= beta.conj() 

138 ptilde += ztilde 

139 else: # First spin 

140 p = z.copy() 

141 ptilde = ztilde.copy() 

142 

143 q = matvec(p) 

144 qtilde = rmatvec(ptilde) 

145 rv = dotprod(ptilde, q) 

146 

147 if rv == 0: 

148 return postprocess(x), -11 

149 

150 alpha = rho_cur / rv 

151 x += alpha*p 

152 r -= alpha*q 

153 rtilde -= alpha.conj()*qtilde 

154 rho_prev = rho_cur 

155 

156 if callback: 

157 callback(x) 

158 

159 else: # for loop exhausted 

160 # Return incomplete progress 

161 return postprocess(x), maxiter 

162 

163 

164@_deprecate_positional_args(version="1.14") 

165def bicgstab(A, b, *, x0=None, tol=_NoValue, maxiter=None, M=None, 

166 callback=None, atol=0., rtol=1e-5): 

167 """Use BIConjugate Gradient STABilized iteration to solve ``Ax = b``. 

168 

169 Parameters 

170 ---------- 

171 A : {sparse matrix, ndarray, LinearOperator} 

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

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

174 produce ``Ax`` and ``A^T x`` using, e.g., 

175 ``scipy.sparse.linalg.LinearOperator``. 

176 b : ndarray 

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

178 x0 : ndarray 

179 Starting guess for the solution. 

180 rtol, atol : float, optional 

181 Parameters for the convergence test. For convergence, 

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

183 The default is ``atol=0.`` and ``rtol=1e-5``. 

184 maxiter : integer 

185 Maximum number of iterations. Iteration will stop after maxiter 

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

187 M : {sparse matrix, ndarray, LinearOperator} 

188 Preconditioner for A. The preconditioner should approximate the 

189 inverse of A. Effective preconditioning dramatically improves the 

190 rate of convergence, which implies that fewer iterations are needed 

191 to reach a given error tolerance. 

192 callback : function 

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

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

195 tol : float, optional, deprecated 

196 

197 .. deprecated 1.12.0 

198 `bicgstab` keyword argument `tol` is deprecated in favor of `rtol` 

199 and will be removed in SciPy 1.14.0. 

200 

201 Returns 

202 ------- 

203 x : ndarray 

204 The converged solution. 

205 info : integer 

206 Provides convergence information: 

207 0 : successful exit 

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

209 <0 : parameter breakdown 

210 

211 Examples 

212 -------- 

213 >>> import numpy as np 

214 >>> from scipy.sparse import csc_matrix 

215 >>> from scipy.sparse.linalg import bicgstab 

216 >>> R = np.array([[4, 2, 0, 1], 

217 ... [3, 0, 0, 2], 

218 ... [0, 1, 1, 1], 

219 ... [0, 2, 1, 0]]) 

220 >>> A = csc_matrix(R) 

221 >>> b = np.array([-1, -0.5, -1, 2]) 

222 >>> x, exit_code = bicgstab(A, b, atol=1e-5) 

223 >>> print(exit_code) # 0 indicates successful convergence 

224 0 

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

226 True 

227 

228 """ 

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

230 bnrm2 = np.linalg.norm(b) 

231 

232 if bnrm2 == 0: 

233 return postprocess(b), 0 

234 

235 atol = _get_atol('bicgstab', b, tol, atol, rtol) 

236 

237 n = len(b) 

238 

239 dotprod = np.vdot if np.iscomplexobj(x) else np.dot 

240 

241 if maxiter is None: 

242 maxiter = n*10 

243 

244 matvec = A.matvec 

245 psolve = M.matvec 

246 

247 # These values make no sense but coming from original Fortran code 

248 # sqrt might have been meant instead. 

249 rhotol = np.finfo(x.dtype.char).eps**2 

250 omegatol = rhotol 

251 

252 # Dummy values to initialize vars, silence linter warnings 

253 rho_prev, omega, alpha, p, v = None, None, None, None, None 

254 

255 r = b - matvec(x) if x.any() else b.copy() 

256 rtilde = r.copy() 

257 

258 for iteration in range(maxiter): 

259 if np.linalg.norm(r) < atol: # Are we done? 

260 return postprocess(x), 0 

261 

262 rho = dotprod(rtilde, r) 

263 if np.abs(rho) < rhotol: # rho breakdown 

264 return postprocess(x), -10 

265 

266 if iteration > 0: 

267 if np.abs(omega) < omegatol: # omega breakdown 

268 return postprocess(x), -11 

269 

270 beta = (rho / rho_prev) * (alpha / omega) 

271 p -= omega*v 

272 p *= beta 

273 p += r 

274 else: # First spin 

275 s = np.empty_like(r) 

276 p = r.copy() 

277 

278 phat = psolve(p) 

279 v = matvec(phat) 

280 rv = dotprod(rtilde, v) 

281 if rv == 0: 

282 return postprocess(x), -11 

283 alpha = rho / rv 

284 r -= alpha*v 

285 s[:] = r[:] 

286 

287 if np.linalg.norm(s) < atol: 

288 x += alpha*phat 

289 return postprocess(x), 0 

290 

291 shat = psolve(s) 

292 t = matvec(shat) 

293 omega = dotprod(t, s) / dotprod(t, t) 

294 x += alpha*phat 

295 x += omega*shat 

296 r -= omega*t 

297 rho_prev = rho 

298 

299 if callback: 

300 callback(x) 

301 

302 else: # for loop exhausted 

303 # Return incomplete progress 

304 return postprocess(x), maxiter 

305 

306 

307@_deprecate_positional_args(version="1.14") 

308def cg(A, b, x0=None, *, tol=_NoValue, maxiter=None, M=None, callback=None, 

309 atol=0., rtol=1e-5): 

310 """Use Conjugate Gradient iteration to solve ``Ax = b``. 

311 

312 Parameters 

313 ---------- 

314 A : {sparse matrix, ndarray, LinearOperator} 

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

316 ``A`` must represent a hermitian, positive definite matrix. 

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

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

319 ``scipy.sparse.linalg.LinearOperator``. 

320 b : ndarray 

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

322 x0 : ndarray 

323 Starting guess for the solution. 

324 rtol, atol : float, optional 

325 Parameters for the convergence test. For convergence, 

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

327 The default is ``atol=0.`` and ``rtol=1e-5``. 

328 maxiter : integer 

329 Maximum number of iterations. Iteration will stop after maxiter 

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

331 M : {sparse matrix, ndarray, LinearOperator} 

332 Preconditioner for A. The preconditioner should approximate the 

333 inverse of A. Effective preconditioning dramatically improves the 

334 rate of convergence, which implies that fewer iterations are needed 

335 to reach a given error tolerance. 

336 callback : function 

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

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

339 tol : float, optional, deprecated 

340 

341 .. deprecated 1.12.0 

342 `cg` keyword argument `tol` is deprecated in favor of `rtol` and 

343 will be removed in SciPy 1.14.0. 

344 

345 Returns 

346 ------- 

347 x : ndarray 

348 The converged solution. 

349 info : integer 

350 Provides convergence information: 

351 0 : successful exit 

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

353 

354 Examples 

355 -------- 

356 >>> import numpy as np 

357 >>> from scipy.sparse import csc_matrix 

358 >>> from scipy.sparse.linalg import cg 

359 >>> P = np.array([[4, 0, 1, 0], 

360 ... [0, 5, 0, 0], 

361 ... [1, 0, 3, 2], 

362 ... [0, 0, 2, 4]]) 

363 >>> A = csc_matrix(P) 

364 >>> b = np.array([-1, -0.5, -1, 2]) 

365 >>> x, exit_code = cg(A, b, atol=1e-5) 

366 >>> print(exit_code) # 0 indicates successful convergence 

367 0 

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

369 True 

370 

371 """ 

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

373 bnrm2 = np.linalg.norm(b) 

374 

375 if bnrm2 == 0: 

376 return postprocess(b), 0 

377 

378 atol = _get_atol('cg', b, tol, atol, rtol) 

379 

380 n = len(b) 

381 

382 if maxiter is None: 

383 maxiter = n*10 

384 

385 dotprod = np.vdot if np.iscomplexobj(x) else np.dot 

386 

387 matvec = A.matvec 

388 psolve = M.matvec 

389 r = b - matvec(x) if x.any() else b.copy() 

390 

391 # Dummy value to initialize var, silences warnings 

392 rho_prev, p = None, None 

393 

394 for iteration in range(maxiter): 

395 if np.linalg.norm(r) < atol: # Are we done? 

396 return postprocess(x), 0 

397 

398 z = psolve(r) 

399 rho_cur = dotprod(r, z) 

400 if iteration > 0: 

401 beta = rho_cur / rho_prev 

402 p *= beta 

403 p += z 

404 else: # First spin 

405 p = np.empty_like(r) 

406 p[:] = z[:] 

407 

408 q = matvec(p) 

409 alpha = rho_cur / dotprod(p, q) 

410 x += alpha*p 

411 r -= alpha*q 

412 rho_prev = rho_cur 

413 

414 if callback: 

415 callback(x) 

416 

417 else: # for loop exhausted 

418 # Return incomplete progress 

419 return postprocess(x), maxiter 

420 

421 

422@_deprecate_positional_args(version="1.14") 

423def cgs(A, b, x0=None, *, tol=_NoValue, maxiter=None, M=None, callback=None, 

424 atol=0., rtol=1e-5): 

425 """Use Conjugate Gradient Squared iteration to solve ``Ax = b``. 

426 

427 Parameters 

428 ---------- 

429 A : {sparse matrix, ndarray, LinearOperator} 

430 The real-valued N-by-N matrix of the linear system. 

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

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

433 ``scipy.sparse.linalg.LinearOperator``. 

434 b : ndarray 

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

436 x0 : ndarray 

437 Starting guess for the solution. 

438 rtol, atol : float, optional 

439 Parameters for the convergence test. For convergence, 

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

441 The default is ``atol=0.`` and ``rtol=1e-5``. 

442 maxiter : integer 

443 Maximum number of iterations. Iteration will stop after maxiter 

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

445 M : {sparse matrix, ndarray, LinearOperator} 

446 Preconditioner for A. The preconditioner should approximate the 

447 inverse of A. Effective preconditioning dramatically improves the 

448 rate of convergence, which implies that fewer iterations are needed 

449 to reach a given error tolerance. 

450 callback : function 

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

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

453 tol : float, optional, deprecated 

454 

455 .. deprecated 1.12.0 

456 `cgs` keyword argument `tol` is deprecated in favor of `rtol` and 

457 will be removed in SciPy 1.14.0. 

458 

459 Returns 

460 ------- 

461 x : ndarray 

462 The converged solution. 

463 info : integer 

464 Provides convergence information: 

465 0 : successful exit 

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

467 <0 : parameter breakdown 

468 

469 Examples 

470 -------- 

471 >>> import numpy as np 

472 >>> from scipy.sparse import csc_matrix 

473 >>> from scipy.sparse.linalg import cgs 

474 >>> R = np.array([[4, 2, 0, 1], 

475 ... [3, 0, 0, 2], 

476 ... [0, 1, 1, 1], 

477 ... [0, 2, 1, 0]]) 

478 >>> A = csc_matrix(R) 

479 >>> b = np.array([-1, -0.5, -1, 2]) 

480 >>> x, exit_code = cgs(A, b) 

481 >>> print(exit_code) # 0 indicates successful convergence 

482 0 

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

484 True 

485 

486 """ 

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

488 bnrm2 = np.linalg.norm(b) 

489 

490 if bnrm2 == 0: 

491 return postprocess(b), 0 

492 

493 atol = _get_atol('cgs', b, tol, atol, rtol) 

494 

495 n = len(b) 

496 

497 dotprod = np.vdot if np.iscomplexobj(x) else np.dot 

498 

499 if maxiter is None: 

500 maxiter = n*10 

501 

502 matvec = A.matvec 

503 psolve = M.matvec 

504 

505 rhotol = np.finfo(x.dtype.char).eps**2 

506 

507 r = b - matvec(x) if x.any() else b.copy() 

508 

509 rtilde = r.copy() 

510 bnorm = np.linalg.norm(b) 

511 if bnorm == 0: 

512 bnorm = 1 

513 

514 # Dummy values to initialize vars, silence linter warnings 

515 rho_prev, p, u, q = None, None, None, None 

516 

517 for iteration in range(maxiter): 

518 rnorm = np.linalg.norm(r) 

519 if rnorm < atol: # Are we done? 

520 return postprocess(x), 0 

521 

522 rho_cur = dotprod(rtilde, r) 

523 if np.abs(rho_cur) < rhotol: # Breakdown case 

524 return postprocess, -10 

525 

526 if iteration > 0: 

527 beta = rho_cur / rho_prev 

528 

529 # u = r + beta * q 

530 # p = u + beta * (q + beta * p); 

531 u[:] = r[:] 

532 u += beta*q 

533 

534 p *= beta 

535 p += q 

536 p *= beta 

537 p += u 

538 

539 else: # First spin 

540 p = r.copy() 

541 u = r.copy() 

542 q = np.empty_like(r) 

543 

544 phat = psolve(p) 

545 vhat = matvec(phat) 

546 rv = dotprod(rtilde, vhat) 

547 

548 if rv == 0: # Dot product breakdown 

549 return postprocess(x), -11 

550 

551 alpha = rho_cur / rv 

552 q[:] = u[:] 

553 q -= alpha*vhat 

554 uhat = psolve(u + q) 

555 x += alpha*uhat 

556 

557 # Due to numerical error build-up the actual residual is computed 

558 # instead of the following two lines that were in the original 

559 # FORTRAN templates, still using a single matvec. 

560 

561 # qhat = matvec(uhat) 

562 # r -= alpha*qhat 

563 r = b - matvec(x) 

564 

565 rho_prev = rho_cur 

566 

567 if callback: 

568 callback(x) 

569 

570 else: # for loop exhausted 

571 # Return incomplete progress 

572 return postprocess(x), maxiter 

573 

574 

575@_deprecate_positional_args(version="1.14") 

576def gmres(A, b, x0=None, *, tol=_NoValue, restart=None, maxiter=None, M=None, 

577 callback=None, restrt=_NoValue, atol=0., callback_type=None, 

578 rtol=1e-5): 

579 """ 

580 Use Generalized Minimal RESidual iteration to solve ``Ax = b``. 

581 

582 Parameters 

583 ---------- 

584 A : {sparse matrix, ndarray, LinearOperator} 

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

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

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

588 ``scipy.sparse.linalg.LinearOperator``. 

589 b : ndarray 

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

591 x0 : ndarray 

592 Starting guess for the solution (a vector of zeros by default). 

593 atol, rtol : float 

594 Parameters for the convergence test. For convergence, 

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

596 The default is ``atol=0.`` and ``rtol=1e-5``. 

597 restart : int, optional 

598 Number of iterations between restarts. Larger values increase 

599 iteration cost, but may be necessary for convergence. 

600 If omitted, ``min(20, n)`` is used. 

601 maxiter : int, optional 

602 Maximum number of iterations (restart cycles). Iteration will stop 

603 after maxiter steps even if the specified tolerance has not been 

604 achieved. See `callback_type`. 

605 M : {sparse matrix, ndarray, LinearOperator} 

606 Inverse of the preconditioner of A. M should approximate the 

607 inverse of A and be easy to solve for (see Notes). Effective 

608 preconditioning dramatically improves the rate of convergence, 

609 which implies that fewer iterations are needed to reach a given 

610 error tolerance. By default, no preconditioner is used. 

611 In this implementation, left preconditioning is used, 

612 and the preconditioned residual is minimized. However, the final 

613 convergence is tested with respect to the ``b - A @ x`` residual. 

614 callback : function 

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

616 as `callback(args)`, where `args` are selected by `callback_type`. 

617 callback_type : {'x', 'pr_norm', 'legacy'}, optional 

618 Callback function argument requested: 

619 - ``x``: current iterate (ndarray), called on every restart 

620 - ``pr_norm``: relative (preconditioned) residual norm (float), 

621 called on every inner iteration 

622 - ``legacy`` (default): same as ``pr_norm``, but also changes the 

623 meaning of `maxiter` to count inner iterations instead of restart 

624 cycles. 

625 

626 This keyword has no effect if `callback` is not set. 

627 restrt : int, optional, deprecated 

628 

629 .. deprecated:: 0.11.0 

630 `gmres` keyword argument `restrt` is deprecated in favor of 

631 `restart` and will be removed in SciPy 1.14.0. 

632 tol : float, optional, deprecated 

633 

634 .. deprecated 1.12.0 

635 `gmres` keyword argument `tol` is deprecated in favor of `rtol` and 

636 will be removed in SciPy 1.14.0 

637 

638 Returns 

639 ------- 

640 x : ndarray 

641 The converged solution. 

642 info : int 

643 Provides convergence information: 

644 0 : successful exit 

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

646 

647 See Also 

648 -------- 

649 LinearOperator 

650 

651 Notes 

652 ----- 

653 A preconditioner, P, is chosen such that P is close to A but easy to solve 

654 for. The preconditioner parameter required by this routine is 

655 ``M = P^-1``. The inverse should preferably not be calculated 

656 explicitly. Rather, use the following template to produce M:: 

657 

658 # Construct a linear operator that computes P^-1 @ x. 

659 import scipy.sparse.linalg as spla 

660 M_x = lambda x: spla.spsolve(P, x) 

661 M = spla.LinearOperator((n, n), M_x) 

662 

663 Examples 

664 -------- 

665 >>> import numpy as np 

666 >>> from scipy.sparse import csc_matrix 

667 >>> from scipy.sparse.linalg import gmres 

668 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float) 

669 >>> b = np.array([2, 4, -1], dtype=float) 

670 >>> x, exitCode = gmres(A, b, atol=1e-5) 

671 >>> print(exitCode) # 0 indicates successful convergence 

672 0 

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

674 True 

675 """ 

676 

677 # Handle the deprecation frenzy 

678 if restrt not in (None, _NoValue) and restart: 

679 raise ValueError("Cannot specify both 'restart' and 'restrt'" 

680 " keywords. Also 'rstrt' is deprecated." 

681 " and will be removed in SciPy 1.14.0. Use " 

682 "'restart' instad.") 

683 if restrt is not _NoValue: 

684 msg = ("'gmres' keyword argument 'restrt' is deprecated " 

685 "in favor of 'restart' and will be removed in SciPy" 

686 " 1.14.0. Until then, if set, 'rstrt' will override 'restart'." 

687 ) 

688 warnings.warn(msg, DeprecationWarning, stacklevel=3) 

689 restart = restrt 

690 

691 if callback is not None and callback_type is None: 

692 # Warn about 'callback_type' semantic changes. 

693 # Probably should be removed only in far future, Scipy 2.0 or so. 

694 msg = ("scipy.sparse.linalg.gmres called without specifying " 

695 "`callback_type`. The default value will be changed in" 

696 " a future release. For compatibility, specify a value " 

697 "for `callback_type` explicitly, e.g., " 

698 "``gmres(..., callback_type='pr_norm')``, or to retain the " 

699 "old behavior ``gmres(..., callback_type='legacy')``" 

700 ) 

701 warnings.warn(msg, category=DeprecationWarning, stacklevel=3) 

702 

703 if callback_type is None: 

704 callback_type = 'legacy' 

705 

706 if callback_type not in ('x', 'pr_norm', 'legacy'): 

707 raise ValueError(f"Unknown callback_type: {callback_type!r}") 

708 

709 if callback is None: 

710 callback_type = None 

711 

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

713 matvec = A.matvec 

714 psolve = M.matvec 

715 n = len(b) 

716 bnrm2 = np.linalg.norm(b) 

717 

718 if bnrm2 == 0: 

719 return postprocess(b), 0 

720 

721 eps = np.finfo(x.dtype.char).eps 

722 

723 dotprod = np.vdot if np.iscomplexobj(x) else np.dot 

724 

725 if maxiter is None: 

726 maxiter = n*10 

727 

728 if restart is None: 

729 restart = 20 

730 restart = min(restart, n) 

731 

732 atol = _get_atol('gmres', b, tol, atol, rtol) 

733 

734 Mb_nrm2 = np.linalg.norm(psolve(b)) 

735 

736 # ==================================================== 

737 # =========== Tolerance control from gh-8400 ========= 

738 # ==================================================== 

739 # Tolerance passed to GMRESREVCOM applies to the inner 

740 # iteration and deals with the left-preconditioned 

741 # residual. 

742 ptol_max_factor = 1. 

743 ptol = Mb_nrm2 * min(ptol_max_factor, atol / bnrm2) 

744 presid = 0. 

745 # ==================================================== 

746 lartg = get_lapack_funcs('lartg', dtype=x.dtype) 

747 

748 # allocate internal variables 

749 v = np.empty([restart+1, n], dtype=x.dtype) 

750 h = np.zeros([restart, restart+1], dtype=x.dtype) 

751 givens = np.zeros([restart, 2], dtype=x.dtype) 

752 

753 # legacy iteration count 

754 inner_iter = 0 

755 

756 for iteration in range(maxiter): 

757 if iteration == 0: 

758 r = b - matvec(x) if x.any() else b.copy() 

759 

760 v[0, :] = psolve(r) 

761 tmp = np.linalg.norm(v[0, :]) 

762 v[0, :] *= (1 / tmp) 

763 # RHS of the Hessenberg problem 

764 S = np.zeros(restart+1, dtype=x.dtype) 

765 S[0] = tmp 

766 

767 breakdown = False 

768 for col in range(restart): 

769 av = matvec(v[col, :]) 

770 w = psolve(av) 

771 

772 # Modified Gram-Schmidt 

773 h0 = np.linalg.norm(w) 

774 for k in range(col+1): 

775 tmp = dotprod(v[k, :], w) 

776 h[col, k] = tmp 

777 w -= tmp*v[k, :] 

778 

779 h1 = np.linalg.norm(w) 

780 h[col, col + 1] = h1 

781 v[col + 1, :] = w[:] 

782 

783 # Exact solution indicator 

784 if h1 <= eps*h0: 

785 h[col, col + 1] = 0 

786 breakdown = True 

787 else: 

788 v[col + 1, :] *= (1 / h1) 

789 

790 # apply past Givens rotations to current h column 

791 for k in range(col): 

792 c, s = givens[k, 0], givens[k, 1] 

793 n0, n1 = h[col, [k, k+1]] 

794 h[col, [k, k + 1]] = [c*n0 + s*n1, -s.conj()*n0 + c*n1] 

795 

796 # get and apply current rotation to h and S 

797 c, s, mag = lartg(h[col, col], h[col, col+1]) 

798 givens[col, :] = [c, s] 

799 h[col, [col, col+1]] = mag, 0 

800 

801 # S[col+1] component is always 0 

802 tmp = -np.conjugate(s)*S[col] 

803 S[[col, col + 1]] = [c*S[col], tmp] 

804 presid = np.abs(tmp) 

805 inner_iter += 1 

806 

807 if callback_type in ('legacy', 'pr_norm'): 

808 callback(presid / bnrm2) 

809 # Legacy behavior 

810 if callback_type == 'legacy' and inner_iter == maxiter: 

811 break 

812 if presid <= ptol or breakdown: 

813 break 

814 

815 # Solve h(col, col) upper triangular system and allow pseudo-solve 

816 # singular cases as in (but without the f2py copies): 

817 # y = trsv(h[:col+1, :col+1].T, S[:col+1]) 

818 

819 if h[col, col] == 0: 

820 S[col] = 0 

821 

822 y = np.zeros([col+1], dtype=x.dtype) 

823 y[:] = S[:col+1] 

824 for k in range(col, 0, -1): 

825 if y[k] != 0: 

826 y[k] /= h[k, k] 

827 tmp = y[k] 

828 y[:k] -= tmp*h[k, :k] 

829 if y[0] != 0: 

830 y[0] /= h[0, 0] 

831 

832 x += y @ v[:col+1, :] 

833 

834 r = b - matvec(x) 

835 rnorm = np.linalg.norm(r) 

836 

837 # Legacy exit 

838 if callback_type == 'legacy' and inner_iter == maxiter: 

839 return postprocess(x), 0 if rnorm <= atol else maxiter 

840 

841 if callback_type == 'x': 

842 callback(x) 

843 

844 if rnorm <= atol: 

845 break 

846 elif breakdown: 

847 # Reached breakdown (= exact solution), but the external 

848 # tolerance check failed. Bail out with failure. 

849 break 

850 elif presid <= ptol: 

851 # Inner loop passed but outer didn't 

852 ptol_max_factor = max(eps, 0.25 * ptol_max_factor) 

853 else: 

854 ptol_max_factor = min(1.0, 1.5 * ptol_max_factor) 

855 

856 ptol = presid * min(ptol_max_factor, atol / rnorm) 

857 

858 info = 0 if (rnorm <= atol) else maxiter 

859 return postprocess(x), info 

860 

861 

862@_deprecate_positional_args(version="1.14") 

863def qmr(A, b, x0=None, *, tol=_NoValue, maxiter=None, M1=None, M2=None, 

864 callback=None, atol=0., rtol=1e-5): 

865 """Use Quasi-Minimal Residual iteration to solve ``Ax = b``. 

866 

867 Parameters 

868 ---------- 

869 A : {sparse matrix, ndarray, LinearOperator} 

870 The real-valued N-by-N matrix of the linear system. 

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

872 produce ``Ax`` and ``A^T x`` using, e.g., 

873 ``scipy.sparse.linalg.LinearOperator``. 

874 b : ndarray 

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

876 x0 : ndarray 

877 Starting guess for the solution. 

878 atol, rtol : float, optional 

879 Parameters for the convergence test. For convergence, 

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

881 The default is ``atol=0.`` and ``rtol=1e-5``. 

882 maxiter : integer 

883 Maximum number of iterations. Iteration will stop after maxiter 

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

885 M1 : {sparse matrix, ndarray, LinearOperator} 

886 Left preconditioner for A. 

887 M2 : {sparse matrix, ndarray, LinearOperator} 

888 Right preconditioner for A. Used together with the left 

889 preconditioner M1. The matrix M1@A@M2 should have better 

890 conditioned than A alone. 

891 callback : function 

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

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

894 tol : float, optional, deprecated 

895 

896 .. deprecated 1.12.0 

897 `qmr` keyword argument `tol` is deprecated in favor of `rtol` and 

898 will be removed in SciPy 1.14.0. 

899 

900 Returns 

901 ------- 

902 x : ndarray 

903 The converged solution. 

904 info : integer 

905 Provides convergence information: 

906 0 : successful exit 

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

908 <0 : parameter breakdown 

909 

910 See Also 

911 -------- 

912 LinearOperator 

913 

914 Examples 

915 -------- 

916 >>> import numpy as np 

917 >>> from scipy.sparse import csc_matrix 

918 >>> from scipy.sparse.linalg import qmr 

919 >>> A = csc_matrix([[3., 2., 0.], [1., -1., 0.], [0., 5., 1.]]) 

920 >>> b = np.array([2., 4., -1.]) 

921 >>> x, exitCode = qmr(A, b, atol=1e-5) 

922 >>> print(exitCode) # 0 indicates successful convergence 

923 0 

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

925 True 

926 """ 

927 A_ = A 

928 A, M, x, b, postprocess = make_system(A, None, x0, b) 

929 bnrm2 = np.linalg.norm(b) 

930 

931 if bnrm2 == 0: