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 817