Hide keyboard shortcuts

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"""Iterative methods for solving linear systems""" 

2 

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

4 

5import warnings 

6import numpy as np 

7 

8from . import _iterative 

9 

10from scipy.sparse.linalg.interface import LinearOperator 

11from .utils import make_system 

12from scipy._lib._util import _aligned_zeros 

13from scipy._lib._threadsafety import non_reentrant 

14 

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

16 

17 

18# Part of the docstring common to all iterative solvers 

19common_doc1 = \ 

20""" 

21Parameters 

22---------- 

23A : {sparse matrix, dense matrix, LinearOperator}""" 

24 

25common_doc2 = \ 

26"""b : {array, matrix} 

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

28 

29Returns 

30------- 

31x : {array, matrix} 

32 The converged solution. 

33info : integer 

34 Provides convergence information: 

35 0 : successful exit 

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

37 <0 : illegal input or breakdown 

38 

39Other Parameters 

40---------------- 

41x0 : {array, matrix} 

42 Starting guess for the solution. 

43tol, atol : float, optional 

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

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

46 a different legacy behavior. 

47 

48 .. warning:: 

49 

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

51 For future compatibility, specify `atol` explicitly. 

52maxiter : integer 

53 Maximum number of iterations. Iteration will stop after maxiter 

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

55M : {sparse matrix, dense matrix, LinearOperator} 

56 Preconditioner for A. The preconditioner should approximate the 

57 inverse of A. Effective preconditioning dramatically improves the 

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

59 to reach a given error tolerance. 

60callback : function 

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

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

63 

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, 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  

135 Examples 

136 -------- 

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 get_residual = lambda: np.linalg.norm(matvec(x) - b) 

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

164 if atol == 'exit': 

165 return postprocess(x), 0 

166 

167 resid = atol 

168 ndx1 = 1 

169 ndx2 = -1 

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

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

172 ijob = 1 

173 info = 0 

174 ftflag = True 

175 iter_ = maxiter 

176 while True: 

177 olditer = iter_ 

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

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

180 if callback is not None and iter_ > olditer: 

181 callback(x) 

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

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

184 if (ijob == -1): 

185 if callback is not None: 

186 callback(x) 

187 break 

188 elif (ijob == 1): 

189 work[slice2] *= sclr2 

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

191 elif (ijob == 2): 

192 work[slice2] *= sclr2 

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

194 elif (ijob == 3): 

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

196 elif (ijob == 4): 

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

198 elif (ijob == 5): 

199 work[slice2] *= sclr2 

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

201 elif (ijob == 6): 

202 if ftflag: 

203 info = -1 

204 ftflag = False 

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

206 ijob = 2 

207 

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

209 # info isn't set appropriately otherwise 

210 info = iter_ 

211 

212 return postprocess(x), info 

213 

214 

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

216 '``Ax = b``.', 

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

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

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

220 '``scipy.sparse.linalg.LinearOperator``.') 

221@non_reentrant() 

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

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

224 

225 n = len(b) 

226 if maxiter is None: 

227 maxiter = n*10 

228 

229 matvec = A.matvec 

230 psolve = M.matvec 

231 ltr = _type_conv[x.dtype.char] 

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

233 

234 get_residual = lambda: np.linalg.norm(matvec(x) - b) 

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

236 if atol == 'exit': 

237 return postprocess(x), 0 

238 

239 resid = atol 

240 ndx1 = 1 

241 ndx2 = -1 

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

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

244 ijob = 1 

245 info = 0 

246 ftflag = True 

247 iter_ = maxiter 

248 while True: 

249 olditer = iter_ 

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

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

252 if callback is not None and iter_ > olditer: 

253 callback(x) 

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

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

256 if (ijob == -1): 

257 if callback is not None: 

258 callback(x) 

259 break 

260 elif (ijob == 1): 

261 work[slice2] *= sclr2 

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

263 elif (ijob == 2): 

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

265 elif (ijob == 3): 

266 work[slice2] *= sclr2 

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

268 elif (ijob == 4): 

269 if ftflag: 

270 info = -1 

271 ftflag = False 

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

273 ijob = 2 

274 

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

276 # info isn't set appropriately otherwise 

277 info = iter_ 

278 

279 return postprocess(x), info 

280 

281 

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

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

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

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

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

287 '``scipy.sparse.linalg.LinearOperator``.') 

288@non_reentrant() 

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

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

291 

292 n = len(b) 

293 if maxiter is None: 

294 maxiter = n*10 

295 

296 matvec = A.matvec 

297 psolve = M.matvec 

298 ltr = _type_conv[x.dtype.char] 

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

300 

301 get_residual = lambda: np.linalg.norm(matvec(x) - b) 

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

303 if atol == 'exit': 

304 return postprocess(x), 0 

305 

306 resid = atol 

307 ndx1 = 1 

308 ndx2 = -1 

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

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

311 ijob = 1 

312 info = 0 

313 ftflag = True 

314 iter_ = maxiter 

315 while True: 

316 olditer = iter_ 

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

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

319 if callback is not None and iter_ > olditer: 

320 callback(x) 

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

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

323 if (ijob == -1): 

324 if callback is not None: 

325 callback(x) 

326 break 

327 elif (ijob == 1): 

328 work[slice2] *= sclr2 

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

330 elif (ijob == 2): 

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

332 elif (ijob == 3): 

333 work[slice2] *= sclr2 

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

335 elif (ijob == 4): 

336 if ftflag: 

337 info = -1 

338 ftflag = False 

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

340 if info == 1 and iter_ > 1: 

341 # recompute residual and recheck, to avoid 

342 # accumulating rounding error 

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

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

345 ijob = 2 

346 

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

348 # info isn't set appropriately otherwise 

349 info = iter_ 

350 

351 return postprocess(x), info 

352 

353 

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

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

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

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

358 '``scipy.sparse.linalg.LinearOperator``.') 

359@non_reentrant() 

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

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

362 

363 n = len(b) 

364 if maxiter is None: 

365 maxiter = n*10 

366 

367 matvec = A.matvec 

368 psolve = M.matvec 

369 ltr = _type_conv[x.dtype.char] 

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

371 

372 get_residual = lambda: np.linalg.norm(matvec(x) - b) 

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

374 if atol == 'exit': 

375 return postprocess(x), 0 

376 

377 resid = atol 

378 ndx1 = 1 

379 ndx2 = -1 

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

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

382 ijob = 1 

383 info = 0 

384 ftflag = True 

385 iter_ = maxiter 

386 while True: 

387 olditer = iter_ 

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

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

390 if callback is not None and iter_ > olditer: 

391 callback(x) 

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

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

394 if (ijob == -1): 

395 if callback is not None: 

396 callback(x) 

397 break 

398 elif (ijob == 1): 

399 work[slice2] *= sclr2 

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

401 elif (ijob == 2): 

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

403 elif (ijob == 3): 

404 work[slice2] *= sclr2 

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

406 elif (ijob == 4): 

407 if ftflag: 

408 info = -1 

409 ftflag = False 

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

411 if info == 1 and iter_ > 1: 

412 # recompute residual and recheck, to avoid 

413 # accumulating rounding error 

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

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

416 ijob = 2 

417 

418 if info == -10: 

419 # termination due to breakdown: check for convergence 

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

421 if ok: 

422 info = 0 

423 

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

425 # info isn't set appropriately otherwise 

426 info = iter_ 

427 

428 return postprocess(x), info 

429 

430 

431@non_reentrant() 

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

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

434 """ 

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

436 

437 Parameters 

438 ---------- 

439 A : {sparse matrix, dense matrix, LinearOperator} 

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

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

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

443 ``scipy.sparse.linalg.LinearOperator``. 

444 b : {array, matrix} 

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

446 

447 Returns 

448 ------- 

449 x : {array, matrix} 

450 The converged solution. 

451 info : int 

452 Provides convergence information: 

453 * 0 : successful exit 

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

455 * <0 : illegal input or breakdown 

456 

457 Other parameters 

458 ---------------- 

459 x0 : {array, matrix} 

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

461 tol, atol : float, optional 

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

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

464 a different legacy behavior. 

465 

466 .. warning:: 

467 

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

469 For future compatibility, specify `atol` explicitly. 

470 restart : int, optional 

471 Number of iterations between restarts. Larger values increase 

472 iteration cost, but may be necessary for convergence. 

473 Default is 20. 

474 maxiter : int, optional 

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

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

477 achieved. 

478 M : {sparse matrix, dense matrix, LinearOperator} 

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

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

481 preconditioning dramatically improves the rate of convergence, 

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

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

484 callback : function 

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

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

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

488 Callback function argument requested: 

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

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

491 called on every inner iteration 

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

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

494 cycles. 

495 restrt : int, optional 

496 DEPRECATED - use `restart` instead. 

497 

498 See Also 

499 -------- 

500 LinearOperator 

501 

502 Notes 

503 ----- 

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

505 for. The preconditioner parameter required by this routine is 

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

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

508 

509 # Construct a linear operator that computes P^-1 * x. 

510 import scipy.sparse.linalg as spla 

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

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

513 

514 Examples 

515 -------- 

516 >>> from scipy.sparse import csc_matrix 

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

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

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

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

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

522 0 

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

524 True 

525 """ 

526 

527 # Change 'restrt' keyword to 'restart' 

528 if restrt is None: 

529 restrt = restart 

530 elif restart is not None: 

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

532 "Preferably use 'restart' only.") 

533 

534 if callback is not None and callback_type is None: 

535 # Warn about 'callback_type' semantic changes. 

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

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

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

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

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

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

542 category=DeprecationWarning, stacklevel=3) 

543 

544 if callback_type is None: 

545 callback_type = 'legacy' 

546 

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

548 raise ValueError("Unknown callback_type: {!r}".format(callback_type)) 

549 

550 if callback is None: 

551 callback_type = 'none' 

552 

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

554 

555 n = len(b) 

556 if maxiter is None: 

557 maxiter = n*10 

558 

559 if restrt is None: 

560 restrt = 20 

561 restrt = min(restrt, n) 

562 

563 matvec = A.matvec 

564 psolve = M.matvec 

565 ltr = _type_conv[x.dtype.char] 

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

567 

568 bnrm2 = np.linalg.norm(b) 

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

570 get_residual = lambda: np.linalg.norm(matvec(x) - b) 

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

572 if atol == 'exit': 

573 return postprocess(x), 0 

574 

575 if bnrm2 == 0: 

576 return postprocess(b), 0 

577 

578 # Tolerance passed to GMRESREVCOM applies to the inner iteration 

579 # and deals with the left-preconditioned residual. 

580 ptol_max_factor = 1.0 

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

582 resid = np.nan 

583 presid = np.nan 

584 ndx1 = 1 

585 ndx2 = -1 

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

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

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

589 ijob = 1 

590 info = 0 

591 ftflag = True 

592 iter_ = maxiter 

593 old_ijob = ijob 

594 first_pass = True 

595 resid_ready = False 

596 iter_num = 1 

597 while True: 

598 olditer = iter_ 

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

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

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

602 callback(x) 

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

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

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

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

607 if resid_ready: 

608 callback(presid / bnrm2) 

609 elif callback_type == 'x': 

610 callback(x) 

611 break 

612 elif (ijob == 1): 

613 work[slice2] *= sclr2 

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

615 elif (ijob == 2): 

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

617 if not first_pass and old_ijob == 3: 

618 resid_ready = True 

619 

620 first_pass = False 

621 elif (ijob == 3): 

622 work[slice2] *= sclr2 

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

624 if resid_ready: 

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

626 callback(presid / bnrm2) 

627 resid_ready = False 

628 iter_num = iter_num+1 

629 

630 elif (ijob == 4): 

631 if ftflag: 

632 info = -1 

633 ftflag = False 

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

635 

636 # Inner loop tolerance control 

637 if info or presid > ptol: 

638 ptol_max_factor = min(1.0, 1.5 * ptol_max_factor) 

639 else: 

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

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

642 

643 if resid != 0: 

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

645 else: 

646 ptol = presid * ptol_max_factor 

647 

648 old_ijob = ijob 

649 ijob = 2 

650 

651 if callback_type == 'legacy': 

652 # Legacy behavior 

653 if iter_num > maxiter: 

654 info = maxiter 

655 break 

656 

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

658 # info isn't set appropriately otherwise 

659 info = maxiter 

660 

661 return postprocess(x), info 

662 

663 

664@non_reentrant() 

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

666 atol=None): 

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

668 

669 Parameters 

670 ---------- 

671 A : {sparse matrix, dense matrix, LinearOperator} 

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

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

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

675 ``scipy.sparse.linalg.LinearOperator``. 

676 b : {array, matrix} 

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

678 

679 Returns 

680 ------- 

681 x : {array, matrix} 

682 The converged solution. 

683 info : integer 

684 Provides convergence information: 

685 0 : successful exit 

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

687 <0 : illegal input or breakdown 

688 

689 Other Parameters 

690 ---------------- 

691 x0 : {array, matrix} 

692 Starting guess for the solution. 

693 tol, atol : float, optional 

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

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

696 a different legacy behavior. 

697 

698 .. warning:: 

699 

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

701 For future compatibility, specify `atol` explicitly. 

702 maxiter : integer 

703 Maximum number of iterations. Iteration will stop after maxiter 

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

705 M1 : {sparse matrix, dense matrix, LinearOperator} 

706 Left preconditioner for A. 

707 M2 : {sparse matrix, dense matrix, LinearOperator} 

708 Right preconditioner for A. Used together with the left 

709 preconditioner M1. The matrix M1*A*M2 should have better 

710 conditioned than A alone. 

711 callback : function 

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

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

714 

715 See Also 

716 -------- 

717 LinearOperator 

718 

719 Examples 

720 -------- 

721 >>> from scipy.sparse import csc_matrix 

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

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

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

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

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

727 0 

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

729 True 

730 """ 

731 A_ = A 

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

733 

734 if M1 is None and M2 is None: 

735 if hasattr(A_,'psolve'): 

736 def left_psolve(b): 

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

738 

739 def right_psolve(b): 

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

741 

742 def left_rpsolve(b): 

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

744 

745 def right_rpsolve(b): 

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

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

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

749 else: 

750 def id(b): 

751 return b 

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

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

754 

755 n = len(b) 

756 if maxiter is None: 

757 maxiter = n*10 

758 

759 ltr = _type_conv[x.dtype.char] 

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

761 

762 get_residual = lambda: np.linalg.norm(A.matvec(x) - b) 

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

764 if atol == 'exit': 

765 return postprocess(x), 0 

766 

767 resid = atol 

768 ndx1 = 1 

769 ndx2 = -1 

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

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

772 ijob = 1 

773 info = 0 

774 ftflag = True 

775 iter_ = maxiter 

776 while True: 

777 olditer = iter_ 

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

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

780 if callback is not None and iter_ > olditer: 

781 callback(x) 

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

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

784 if (ijob == -1): 

785 if callback is not None: 

786 callback(x) 

787 break 

788 elif (ijob == 1): 

789 work[slice2] *= sclr2 

790 work[slice2] += sclr1*A.matvec(work[slice1]) 

791 elif (ijob == 2): 

792 work[slice2] *= sclr2 

793 work[slice2] += sclr1*A.rmatvec(work[slice1]) 

794 elif (ijob == 3): 

795 work[slice1] = M1.matvec(work[slice2]) 

796 elif (ijob == 4): 

797 work[slice1] = M2.matvec(work[slice2]) 

798 elif (ijob == 5): 

799 work[slice1] = M1.rmatvec(work[slice2]) 

800 elif (ijob == 6): 

801 work[slice1] = M2.rmatvec(work[slice2]) 

802 elif (ijob == 7): 

803 work[slice2] *= sclr2 

804 work[slice2] += sclr1*A.matvec(x) 

805 elif (ijob == 8): 

806 if ftflag: 

807 info = -1 

808 ftflag = False 

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

810 ijob = 2 

811 

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

813 # info isn't set appropriately otherwise 

814 info = iter_ 

815 

816 return postprocess(x), info