Coverage for /usr/lib/python3/dist-packages/scipy/sparse/linalg/_isolve/iterative.py: 8%

1"""Iterative methods for solving linear systems""" 

2 

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

4 

5import warnings 

6from textwrap import dedent 

7import numpy as np 

8 

9from . import _iterative 

10 

11from scipy.sparse.linalg._interface import LinearOperator 

12from .utils import make_system 

13from scipy._lib._util import _aligned_zeros 

14from scipy._lib._threadsafety import non_reentrant 

15 

16_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'} 

17 

18 

19# Part of the docstring common to all iterative solvers 

20common_doc1 = \ 

21""" 

22Parameters 

23---------- 

24A : {sparse matrix, ndarray, LinearOperator}""" 

25 

26common_doc2 = \ 

27"""b : ndarray 

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

29 

30Returns 

31------- 

32x : ndarray 

33 The converged solution. 

34info : integer 

35 Provides convergence information: 

36 0 : successful exit 

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

38 <0 : illegal input or breakdown 

39 

40Other Parameters 

41---------------- 

42x0 : ndarray 

43 Starting guess for the solution. 

44tol, atol : float, optional 

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

46 The default for ``atol`` is ``'legacy'``, which emulates 

47 a different legacy behavior. 

48 

49 .. warning:: 

50 

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

52 For future compatibility, specify `atol` explicitly. 

53maxiter : integer 

54 Maximum number of iterations. Iteration will stop after maxiter 

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

56M : {sparse matrix, ndarray, LinearOperator} 

57 Preconditioner for A. The preconditioner should approximate the 

58 inverse of A. Effective preconditioning dramatically improves the 

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

60 to reach a given error tolerance. 

61callback : function 

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

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

64""" 

65 

66 

67def _stoptest(residual, atol): 

68 """ 

69 Successful termination condition for the solvers. 

70 """ 

71 resid = np.linalg.norm(residual) 

72 if resid <= atol: 

73 return resid, 1 

74 else: 

75 return resid, 0 

76 

77 

78def _get_atol(tol, atol, bnrm2, get_residual, routine_name): 

79 """ 

80 Parse arguments for absolute tolerance in termination condition. 

81 

82 Parameters 

83 ---------- 

84 tol, atol : object 

85 The arguments passed into the solver routine by user. 

86 bnrm2 : float 

87 2-norm of the rhs vector. 

88 get_residual : callable 

89 Callable ``get_residual()`` that returns the initial value of 

90 the residual. 

91 routine_name : str 

92 Name of the routine. 

93 """ 

94 

95 if atol is None: 

96 warnings.warn("scipy.sparse.linalg.{name} called without specifying `atol`. " 

97 "The default value will be changed in a future release. " 

98 "For compatibility, specify a value for `atol` explicitly, e.g., " 

99 "``{name}(..., atol=0)``, or to retain the old behavior " 

100 "``{name}(..., atol='legacy')``".format(name=routine_name), 

101 category=DeprecationWarning, stacklevel=4) 

102 atol = 'legacy' 

103 

104 tol = float(tol) 

105 

106 if atol == 'legacy': 

107 # emulate old legacy behavior 

108 resid = get_residual() 

109 if resid <= tol: 

110 return 'exit' 

111 if bnrm2 == 0: 

112 return tol 

113 else: 

114 return tol * float(bnrm2) 

115 else: 

116 return max(float(atol), tol * float(bnrm2)) 

117 

118 

119def set_docstring(header, Ainfo, footer='', atol_default='0'): 

120 def combine(fn): 

121 fn.__doc__ = '\n'.join((header, common_doc1, 

122 ' ' + Ainfo.replace('\n', '\n '), 

123 common_doc2, dedent(footer))) 

124 return fn 

125 return combine 

126 

127 

128@set_docstring('Use BIConjugate Gradient iteration to solve ``Ax = b``.', 

129 'The real or complex N-by-N matrix of the linear system.\n' 

130 'Alternatively, ``A`` can be a linear operator which can\n' 

131 'produce ``Ax`` and ``A^T x`` using, e.g.,\n' 

132 '``scipy.sparse.linalg.LinearOperator``.', 

133 footer="""\ 

134 Examples 

135 -------- 

136 >>> import numpy as np 

137 >>> from scipy.sparse import csc_matrix 

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

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

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

141 >>> x, exitCode = bicg(A, b) 

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

143 0 

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

145 True 

146 

147 """ 

148 ) 

149@non_reentrant() 

150def bicg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None): 

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

152 

153 n = len(b) 

154 if maxiter is None: 

155 maxiter = n*10 

156 

157 matvec, rmatvec = A.matvec, A.rmatvec 

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

159 ltr = _type_conv[x.dtype.char] 

160 revcom = getattr(_iterative, ltr + 'bicgrevcom') 

161 

162 def get_residual(): 

163 return np.linalg.norm(matvec(x) - b) 

164 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'bicg') 

165 if atol == 'exit': 

166 return postprocess(x), 0 

167 

168 resid = atol 

169 ndx1 = 1 

170 ndx2 = -1 

171 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

172 work = _aligned_zeros(6*n,dtype=x.dtype) 

173 ijob = 1 

174 info = 0 

175 ftflag = True 

176 iter_ = maxiter 

177 while True: 

178 olditer = iter_ 

179 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ 

180 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) 

181 if callback is not None and iter_ > olditer: 

182 callback(x) 

183 slice1 = slice(ndx1-1, ndx1-1+n) 

184 slice2 = slice(ndx2-1, ndx2-1+n) 

185 if (ijob == -1): 

186 if callback is not None: 

187 callback(x) 

188 break 

189 elif (ijob == 1): 

190 work[slice2] *= sclr2 

191 work[slice2] += sclr1*matvec(work[slice1]) 

192 elif (ijob == 2): 

193 work[slice2] *= sclr2 

194 work[slice2] += sclr1*rmatvec(work[slice1]) 

195 elif (ijob == 3): 

196 work[slice1] = psolve(work[slice2]) 

197 elif (ijob == 4): 

198 work[slice1] = rpsolve(work[slice2]) 

199 elif (ijob == 5): 

200 work[slice2] *= sclr2 

201 work[slice2] += sclr1*matvec(x) 

202 elif (ijob == 6): 

203 if ftflag: 

204 info = -1 

205 ftflag = False 

206 resid, info = _stoptest(work[slice1], atol) 

207 ijob = 2 

208 

209 if info > 0 and iter_ == maxiter and not (resid <= atol): 

210 # info isn't set appropriately otherwise 

211 info = iter_ 

212 

213 return postprocess(x), info 

214 

215 

216@set_docstring('Use BIConjugate Gradient STABilized iteration to solve ' 

217 '``Ax = b``.', 

218 'The real or complex N-by-N matrix of the linear system.\n' 

219 'Alternatively, ``A`` can be a linear operator which can\n' 

220 'produce ``Ax`` using, e.g.,\n' 

221 '``scipy.sparse.linalg.LinearOperator``.', 

222 footer="""\ 

223 Examples 

224 -------- 

225 >>> import numpy as np 

226 >>> from scipy.sparse import csc_matrix 

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

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

229 ... [3, 0, 0, 2], 

230 ... [0, 1, 1, 1], 

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

232 >>> A = csc_matrix(R) 

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

234 >>> x, exit_code = bicgstab(A, b) 

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

236 0 

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

238 True 

239 """) 

240@non_reentrant() 

241def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None): 

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

243 

244 n = len(b) 

245 if maxiter is None: 

246 maxiter = n*10 

247 

248 matvec = A.matvec 

249 psolve = M.matvec 

250 ltr = _type_conv[x.dtype.char] 

251 revcom = getattr(_iterative, ltr + 'bicgstabrevcom') 

252 

253 def get_residual(): 

254 return np.linalg.norm(matvec(x) - b) 

255 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'bicgstab') 

256 if atol == 'exit': 

257 return postprocess(x), 0 

258 

259 resid = atol 

260 ndx1 = 1 

261 ndx2 = -1 

262 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

263 work = _aligned_zeros(7*n,dtype=x.dtype) 

264 ijob = 1 

265 info = 0 

266 ftflag = True 

267 iter_ = maxiter 

268 while True: 

269 olditer = iter_ 

270 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ 

271 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) 

272 if callback is not None and iter_ > olditer: 

273 callback(x) 

274 slice1 = slice(ndx1-1, ndx1-1+n) 

275 slice2 = slice(ndx2-1, ndx2-1+n) 

276 if (ijob == -1): 

277 if callback is not None: 

278 callback(x) 

279 break 

280 elif (ijob == 1): 

281 work[slice2] *= sclr2 

282 work[slice2] += sclr1*matvec(work[slice1]) 

283 elif (ijob == 2): 

284 work[slice1] = psolve(work[slice2]) 

285 elif (ijob == 3): 

286 work[slice2] *= sclr2 

287 work[slice2] += sclr1*matvec(x) 

288 elif (ijob == 4): 

289 if ftflag: 

290 info = -1 

291 ftflag = False 

292 resid, info = _stoptest(work[slice1], atol) 

293 ijob = 2 

294 

295 if info > 0 and iter_ == maxiter and not (resid <= atol): 

296 # info isn't set appropriately otherwise 

297 info = iter_ 

298 

299 return postprocess(x), info 

300 

301 

302@set_docstring('Use Conjugate Gradient iteration to solve ``Ax = b``.', 

303 'The real or complex N-by-N matrix of the linear system.\n' 

304 '``A`` must represent a hermitian, positive definite matrix.\n' 

305 'Alternatively, ``A`` can be a linear operator which can\n' 

306 'produce ``Ax`` using, e.g.,\n' 

307 '``scipy.sparse.linalg.LinearOperator``.', 

308 footer="""\ 

309 Examples 

310 -------- 

311 >>> import numpy as np 

312 >>> from scipy.sparse import csc_matrix 

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

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

315 ... [0, 5, 0, 0], 

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

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

318 >>> A = csc_matrix(P) 

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

320 >>> x, exit_code = cg(A, b) 

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

322 0 

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

324 True 

325 

326 """) 

327@non_reentrant() 

328def cg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None): 

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

330 

331 n = len(b) 

332 if maxiter is None: 

333 maxiter = n*10 

334 

335 matvec = A.matvec 

336 psolve = M.matvec 

337 ltr = _type_conv[x.dtype.char] 

338 revcom = getattr(_iterative, ltr + 'cgrevcom') 

339 

340 def get_residual(): 

341 return np.linalg.norm(matvec(x) - b) 

342 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'cg') 

343 if atol == 'exit': 

344 return postprocess(x), 0 

345 

346 resid = atol 

347 ndx1 = 1 

348 ndx2 = -1 

349 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

350 work = _aligned_zeros(4*n,dtype=x.dtype) 

351 ijob = 1 

352 info = 0 

353 ftflag = True 

354 iter_ = maxiter 

355 while True: 

356 olditer = iter_ 

357 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ 

358 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) 

359 if callback is not None and iter_ > olditer: 

360 callback(x) 

361 slice1 = slice(ndx1-1, ndx1-1+n) 

362 slice2 = slice(ndx2-1, ndx2-1+n) 

363 if (ijob == -1): 

364 if callback is not None: 

365 callback(x) 

366 break 

367 elif (ijob == 1): 

368 work[slice2] *= sclr2 

369 work[slice2] += sclr1*matvec(work[slice1]) 

370 elif (ijob == 2): 

371 work[slice1] = psolve(work[slice2]) 

372 elif (ijob == 3): 

373 work[slice2] *= sclr2 

374 work[slice2] += sclr1*matvec(x) 

375 elif (ijob == 4): 

376 if ftflag: 

377 info = -1 

378 ftflag = False 

379 resid, info = _stoptest(work[slice1], atol) 

380 if info == 1 and iter_ > 1: 

381 # recompute residual and recheck, to avoid 

382 # accumulating rounding error 

383 work[slice1] = b - matvec(x) 

384 resid, info = _stoptest(work[slice1], atol) 

385 ijob = 2 

386 

387 if info > 0 and iter_ == maxiter and not (resid <= atol): 

388 # info isn't set appropriately otherwise 

389 info = iter_ 

390 

391 return postprocess(x), info 

392 

393 

394@set_docstring('Use Conjugate Gradient Squared iteration to solve ``Ax = b``.', 

395 'The real-valued N-by-N matrix of the linear system.\n' 

396 'Alternatively, ``A`` can be a linear operator which can\n' 

397 'produce ``Ax`` using, e.g.,\n' 

398 '``scipy.sparse.linalg.LinearOperator``.', 

399 footer="""\ 

400 Examples 

401 -------- 

402 >>> import numpy as np 

403 >>> from scipy.sparse import csc_matrix 

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

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

406 ... [3, 0, 0, 2], 

407 ... [0, 1, 1, 1], 

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

409 >>> A = csc_matrix(R) 

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

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

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

413 0 

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

415 True 

416 """ 

417 ) 

418@non_reentrant() 

419def cgs(A, b, x0=None, tol=1e-5, maxiter=None, M=None, callback=None, atol=None): 

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

421 

422 n = len(b) 

423 if maxiter is None: 

424 maxiter = n*10 

425 

426 matvec = A.matvec 

427 psolve = M.matvec 

428 ltr = _type_conv[x.dtype.char] 

429 revcom = getattr(_iterative, ltr + 'cgsrevcom') 

430 

431 def get_residual(): 

432 return np.linalg.norm(matvec(x) - b) 

433 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'cgs') 

434 if atol == 'exit': 

435 return postprocess(x), 0 

436 

437 resid = atol 

438 ndx1 = 1 

439 ndx2 = -1 

440 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

441 work = _aligned_zeros(7*n,dtype=x.dtype) 

442 ijob = 1 

443 info = 0 

444 ftflag = True 

445 iter_ = maxiter 

446 while True: 

447 olditer = iter_ 

448 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ 

449 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) 

450 if callback is not None and iter_ > olditer: 

451 callback(x) 

452 slice1 = slice(ndx1-1, ndx1-1+n) 

453 slice2 = slice(ndx2-1, ndx2-1+n) 

454 if (ijob == -1): 

455 if callback is not None: 

456 callback(x) 

457 break 

458 elif (ijob == 1): 

459 work[slice2] *= sclr2 

460 work[slice2] += sclr1*matvec(work[slice1]) 

461 elif (ijob == 2): 

462 work[slice1] = psolve(work[slice2]) 

463 elif (ijob == 3): 

464 work[slice2] *= sclr2 

465 work[slice2] += sclr1*matvec(x) 

466 elif (ijob == 4): 

467 if ftflag: 

468 info = -1 

469 ftflag = False 

470 resid, info = _stoptest(work[slice1], atol) 

471 if info == 1 and iter_ > 1: 

472 # recompute residual and recheck, to avoid 

473 # accumulating rounding error 

474 work[slice1] = b - matvec(x) 

475 resid, info = _stoptest(work[slice1], atol) 

476 ijob = 2 

477 

478 if info == -10: 

479 # termination due to breakdown: check for convergence 

480 resid, ok = _stoptest(b - matvec(x), atol) 

481 if ok: 

482 info = 0 

483 

484 if info > 0 and iter_ == maxiter and not (resid <= atol): 

485 # info isn't set appropriately otherwise 

486 info = iter_ 

487 

488 return postprocess(x), info 

489 

490 

491@non_reentrant() 

492def gmres(A, b, x0=None, tol=1e-5, restart=None, maxiter=None, M=None, callback=None, 

493 restrt=None, atol=None, callback_type=None): 

494 """ 

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

496 

497 Parameters 

498 ---------- 

499 A : {sparse matrix, ndarray, LinearOperator} 

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

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

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

503 ``scipy.sparse.linalg.LinearOperator``. 

504 b : ndarray 

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

506 

507 Returns 

508 ------- 

509 x : ndarray 

510 The converged solution. 

511 info : int 

512 Provides convergence information: 

513 * 0 : successful exit 

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

515 * <0 : illegal input or breakdown 

516 

517 Other parameters 

518 ---------------- 

519 x0 : ndarray 

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

521 tol, atol : float, optional 

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

523 The default for ``atol`` is ``'legacy'``, which emulates 

524 a different legacy behavior. 

525 

526 .. warning:: 

527 

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

529 For future compatibility, specify `atol` explicitly. 

530 restart : int, optional 

531 Number of iterations between restarts. Larger values increase 

532 iteration cost, but may be necessary for convergence. 

533 Default is 20. 

534 maxiter : int, optional 

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

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

537 achieved. 

538 M : {sparse matrix, ndarray, LinearOperator} 

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

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

541 preconditioning dramatically improves the rate of convergence, 

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

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

544 In this implementation, left preconditioning is used, 

545 and the preconditioned residual is minimized. 

546 callback : function 

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

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

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

550 Callback function argument requested: 

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

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

553 called on every inner iteration 

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

555 meaning of 'maxiter' to count inner iterations instead of restart 

556 cycles. 

557 restrt : int, optional, deprecated 

558 

559 .. deprecated:: 0.11.0 

560 `gmres` keyword argument `restrt` is deprecated infavour of 

561 `restart` and will be removed in SciPy 1.12.0. 

562 

563 See Also 

564 -------- 

565 LinearOperator 

566 

567 Notes 

568 ----- 

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

570 for. The preconditioner parameter required by this routine is 

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

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

573 

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

575 import scipy.sparse.linalg as spla 

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

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

578 

579 Examples 

580 -------- 

581 >>> import numpy as np 

582 >>> from scipy.sparse import csc_matrix 

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

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

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

586 >>> x, exitCode = gmres(A, b) 

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

588 0 

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

590 True 

591 """ 

592 

593 # Change 'restrt' keyword to 'restart' 

594 if restrt is None: 

595 restrt = restart 

596 elif restart is not None: 

597 raise ValueError("Cannot specify both restart and restrt keywords. " 

598 "Preferably use 'restart' only.") 

599 else: 

600 msg = ("'gmres' keyword argument 'restrt' is deprecated infavour of " 

601 "'restart' and will be removed in SciPy 1.12.0.") 

602 warnings.warn(msg, DeprecationWarning, stacklevel=2) 

603 

604 if callback is not None and callback_type is None: 

605 # Warn about 'callback_type' semantic changes. 

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

607 warnings.warn("scipy.sparse.linalg.gmres called without specifying `callback_type`. " 

608 "The default value will be changed in a future release. " 

609 "For compatibility, specify a value for `callback_type` explicitly, e.g., " 

610 "``{name}(..., callback_type='pr_norm')``, or to retain the old behavior " 

611 "``{name}(..., callback_type='legacy')``", 

612 category=DeprecationWarning, stacklevel=3) 

613 

614 if callback_type is None: 

615 callback_type = 'legacy' 

616 

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

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

619 

620 if callback is None: 

621 callback_type = 'none' 

622 

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

624 

625 n = len(b) 

626 if maxiter is None: 

627 maxiter = n*10 

628 

629 if restrt is None: 

630 restrt = 20 

631 restrt = min(restrt, n) 

632 

633 matvec = A.matvec 

634 psolve = M.matvec 

635 ltr = _type_conv[x.dtype.char] 

636 revcom = getattr(_iterative, ltr + 'gmresrevcom') 

637 

638 bnrm2 = np.linalg.norm(b) 

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

640 def get_residual(): 

641 return np.linalg.norm(matvec(x) - b) 

642 atol = _get_atol(tol, atol, bnrm2, get_residual, 'gmres') 

643 if atol == 'exit': 

644 return postprocess(x), 0 

645 

646 if bnrm2 == 0: 

647 return postprocess(b), 0 

648 

649 # Tolerance passed to GMRESREVCOM applies to the inner iteration 

650 # and deals with the left-preconditioned residual. 

651 ptol_max_factor = 1.0 

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

653 resid = np.nan 

654 presid = np.nan 

655 ndx1 = 1 

656 ndx2 = -1 

657 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

658 work = _aligned_zeros((6+restrt)*n,dtype=x.dtype) 

659 work2 = _aligned_zeros((restrt+1)*(2*restrt+2),dtype=x.dtype) 

660 ijob = 1 

661 info = 0 

662 ftflag = True 

663 iter_ = maxiter 

664 old_ijob = ijob 

665 first_pass = True 

666 resid_ready = False 

667 iter_num = 1 

668 while True: 

669 olditer = iter_ 

670 x, iter_, presid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ 

671 revcom(b, x, restrt, work, work2, iter_, presid, info, ndx1, ndx2, ijob, ptol) 

672 if callback_type == 'x' and iter_ != olditer: 

673 callback(x) 

674 slice1 = slice(ndx1-1, ndx1-1+n) 

675 slice2 = slice(ndx2-1, ndx2-1+n) 

676 if (ijob == -1): # gmres success, update last residual 

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

678 if resid_ready: 

679 callback(presid / bnrm2) 

680 elif callback_type == 'x': 

681 callback(x) 

682 break 

683 elif (ijob == 1): 

684 work[slice2] *= sclr2 

685 work[slice2] += sclr1*matvec(x) 

686 elif (ijob == 2): 

687 work[slice1] = psolve(work[slice2]) 

688 if not first_pass and old_ijob == 3: 

689 resid_ready = True 

690 

691 first_pass = False 

692 elif (ijob == 3): 

693 work[slice2] *= sclr2 

694 work[slice2] += sclr1*matvec(work[slice1]) 

695 if resid_ready: 

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

697 callback(presid / bnrm2) 

698 resid_ready = False 

699 iter_num = iter_num+1 

700 

701 elif (ijob == 4): 

702 if ftflag: 

703 info = -1 

704 ftflag = False 

705 resid, info = _stoptest(work[slice1], atol) 

706 

707 # Inner loop tolerance control 

708 if info or presid > ptol: 

709 ptol_max_factor = min(1.0, 1.5 * ptol_max_factor) 

710 else: 

711 # Inner loop tolerance OK, but outer loop not. 

712 ptol_max_factor = max(1e-16, 0.25 * ptol_max_factor) 

713 

714 if resid != 0: 

715 ptol = presid * min(ptol_max_factor, atol / resid) 

716 else: 

717 ptol = presid * ptol_max_factor 

718 

719 old_ijob = ijob 

720 ijob = 2 

721 

722 if callback_type == 'legacy': 

723 # Legacy behavior 

724 if iter_num > maxiter: 

725 info = maxiter 

726 break 

727 

728 if info >= 0 and not (resid <= atol): 

729 # info isn't set appropriately otherwise 

730 info = maxiter 

731 

732 return postprocess(x), info 

733 

734 

735@non_reentrant() 

736def qmr(A, b, x0=None, tol=1e-5, maxiter=None, M1=None, M2=None, callback=None, 

737 atol=None): 

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

739 

740 Parameters 

741 ---------- 

742 A : {sparse matrix, ndarray, LinearOperator} 

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

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

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

746 ``scipy.sparse.linalg.LinearOperator``. 

747 b : ndarray 

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

749 

750 Returns 

751 ------- 

752 x : ndarray 

753 The converged solution. 

754 info : integer 

755 Provides convergence information: 

756 0 : successful exit 

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

758 <0 : illegal input or breakdown 

759 

760 Other Parameters 

761 ---------------- 

762 x0 : ndarray 

763 Starting guess for the solution. 

764 tol, atol : float, optional 

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

766 The default for ``atol`` is ``'legacy'``, which emulates 

767 a different legacy behavior. 

768 

769 .. warning:: 

770 

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

772 For future compatibility, specify `atol` explicitly. 

773 maxiter : integer 

774 Maximum number of iterations. Iteration will stop after maxiter 

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

776 M1 : {sparse matrix, ndarray, LinearOperator} 

777 Left preconditioner for A. 

778 M2 : {sparse matrix, ndarray, LinearOperator} 

779 Right preconditioner for A. Used together with the left 

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

781 conditioned than A alone. 

782 callback : function 

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

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

785 

786 See Also 

787 -------- 

788 LinearOperator 

789 

790 Examples 

791 -------- 

792 >>> import numpy as np 

793 >>> from scipy.sparse import csc_matrix 

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

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

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

797 >>> x, exitCode = qmr(A, b) 

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

799 0 

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

801 True 

802 """ 

803 A_ = A 

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

805 

806 if M1 is None and M2 is None: 

807 if hasattr(A_,'psolve'): 

808 def left_psolve(b): 

809 return A_.psolve(b,'left') 

810 

811 def right_psolve(b): 

812 return A_.psolve(b,'right') 

813 

814 def left_rpsolve(b): 

815 return A_.rpsolve(b,'left') 

816 

817 def right_rpsolve(b): 

818 return A_.rpsolve(b,'right') 

819 M1 = LinearOperator(A.shape, matvec=left_psolve, rmatvec=left_rpsolve) 

820 M2 = LinearOperator(A.shape, matvec=right_psolve, rmatvec=right_rpsolve) 

821 else: 

822 def id(b): 

823 return b 

824 M1 = LinearOperator(A.shape, matvec=id, rmatvec=id) 

825 M2 = LinearOperator(A.shape, matvec=id, rmatvec=id) 

826 

827 n = len(b) 

828 if maxiter is None: 

829 maxiter = n*10 

830 

831 ltr = _type_conv[x.dtype.char] 

832 revcom = getattr(_iterative, ltr + 'qmrrevcom') 

833 

834 def get_residual(): 

835 return np.linalg.norm(A.matvec(x) - b) 

836 atol = _get_atol(tol, atol, np.linalg.norm(b), get_residual, 'qmr') 

837 if atol == 'exit': 

838 return postprocess(x), 0 

839 

840 resid = atol 

841 ndx1 = 1 

842 ndx2 = -1 

843 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

844 work = _aligned_zeros(11*n,x.dtype) 

845 ijob = 1 

846 info = 0 

847 ftflag = True 

848 iter_ = maxiter 

849 while True: 

850 olditer = iter_ 

851 x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \ 

852 revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob) 

853 if callback is not None and iter_ > olditer: 

854 callback(x) 

855 slice1 = slice(ndx1-1, ndx1-1+n) 

856 slice2 = slice(ndx2-1, ndx2-1+n) 

857 if (ijob == -1): 

858 if callback is not None: 

859 callback(x) 

860 break