code: """
  Copyright (C) 2010 David Fong and Michael Saunders
  
  LSMR uses an iterative method.
  
  07 Jun 2010: Documentation updated
  03 Jun 2010: First release version in Python
  
  David Chin-lung Fong            clfong@stanford.edu
  Institute for Computational and Mathematical Engineering
  Stanford University
  
  Michael Saunders                saunders@stanford.edu
  Systems Optimization Laboratory
  Dept of MS&E, Stanford University.
  
  """
  
  from __future__ import division, print_function, absolute_import
  
  __all__ = ['lsmr']
  
  from numpy import zeros, infty
  from numpy.linalg import norm
  from math import sqrt
  from scipy.sparse.linalg.interface import aslinearoperator
  
  from .lsqr import _sym_ortho
  
  def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
           maxiter=None, show=False):
      """Iterative solver for least-squares problems.
  
      lsmr solves the system of linear equations ``Ax = b``. If the system
      is inconsistent, it solves the least-squares problem ``min ||b - Ax||_2``.
      A is a rectangular matrix of dimension m-by-n, where all cases are
      allowed: m = n, m > n, or m < n. B is a vector of length m.
      The matrix A may be dense or sparse (usually sparse).
  
      .. versionadded:: 0.11.0
  
      Parameters
      ----------
      A : {matrix, sparse matrix, ndarray, LinearOperator}
          Matrix A in the linear system.
      b : (m,) ndarray
          Vector b in the linear system.
      damp : float
          Damping factor for regularized least-squares. `lsmr` solves
          the regularized least-squares problem::
  
           min ||(b) - (  A   )x||
               ||(0)   (damp*I) ||_2
  
          where damp is a scalar.  If damp is None or 0, the system
          is solved without regularization.
      atol, btol : float
          Stopping tolerances. `lsmr` continues iterations until a
          certain backward error estimate is smaller than some quantity
          depending on atol and btol.  Let ``r = b - Ax`` be the
          residual vector for the current approximate solution ``x``.
          If ``Ax = b`` seems to be consistent, ``lsmr`` terminates
          when ``norm(r) <= atol * norm(A) * norm(x) + btol * norm(b)``.
          Otherwise, lsmr terminates when ``norm(A^{T} r) <=
          atol * norm(A) * norm(r)``.  If both tolerances are 1.0e-6 (say),
          the final ``norm(r)`` should be accurate to about 6
          digits. (The final x will usually have fewer correct digits,
          depending on ``cond(A)`` and the size of LAMBDA.)  If `atol`
          or `btol` is None, a default value of 1.0e-6 will be used.
          Ideally, they should be estimates of the relative error in the
          entries of A and B respectively.  For example, if the entries
          of `A` have 7 correct digits, set atol = 1e-7. This prevents
          the algorithm from doing unnecessary work beyond the
          uncertainty of the input data.
      conlim : float
          `lsmr` terminates if an estimate of ``cond(A)`` exceeds
          `conlim`.  For compatible systems ``Ax = b``, conlim could be
          as large as 1.0e+12 (say).  For least-squares problems,
          `conlim` should be less than 1.0e+8. If `conlim` is None, the
          default value is 1e+8.  Maximum precision can be obtained by
          setting ``atol = btol = conlim = 0``, but the number of
          iterations may then be excessive.
      maxiter : int
          `lsmr` terminates if the number of iterations reaches
          `maxiter`.  The default is ``maxiter = min(m, n)``.  For
          ill-conditioned systems, a larger value of `maxiter` may be
          needed.
      show : bool
          Print iterations logs if ``show=True``.
  
      Returns
      -------
      x : ndarray of float
          Least-square solution returned.
      istop : int
          istop gives the reason for stopping::
  
            istop   = 0 means x=0 is a solution.
                    = 1 means x is an approximate solution to A*x = B,
                        according to atol and btol.
                    = 2 means x approximately solves the least-squares problem
                        according to atol.
                    = 3 means COND(A) seems to be greater than CONLIM.
                    = 4 is the same as 1 with atol = btol = eps (machine
                        precision)
                    = 5 is the same as 2 with atol = eps.
                    = 6 is the same as 3 with CONLIM = 1/eps.
                    = 7 means ITN reached maxiter before the other stopping
                        conditions were satisfied.
  
      itn : int
          Number of iterations used.
      normr : float
          ``norm(b-Ax)``
      normar : float
          ``norm(A^T (b - Ax))``
      norma : float
          ``norm(A)``
      conda : float
          Condition number of A.
      normx : float
          ``norm(x)``
  
      References
      ----------
      .. [1] D. C.-L. Fong and M. A. Saunders,
             "LSMR: An iterative algorithm for sparse least-squares problems",
             SIAM J. Sci. Comput., vol. 33, pp. 2950-2971, 2011.
             http://arxiv.org/abs/1006.0758
      .. [2] LSMR Software, http://www.stanford.edu/~clfong/lsmr.html
  
      """
  
      A = aslinearoperator(A)
      b = b.squeeze()
  
      msg=('The exact solution is  x = 0                              ',
           'Ax - b is small enough, given atol, btol                  ',
           'The least-squares solution is good enough, given atol     ',
           'The estimate of cond(Abar) has exceeded conlim            ',
           'Ax - b is small enough for this machine                   ',
           'The least-squares solution is good enough for this machine',
           'Cond(Abar) seems to be too large for this machine         ',
           'The iteration limit has been reached                      ')
  
      hdg1 = '   itn      x(1)       norm r    norm A''r'
      hdg2 = ' compatible   LS      norm A   cond A'
      pfreq = 20   # print frequency (for repeating the heading)
      pcount = 0   # print counter
  
      m, n = A.shape
  
      # stores the num of singular values
      minDim = min([m, n])
  
      if maxiter is None:
          maxiter = minDim
  
      if show:
          print(' ')
          print('LSMR            Least-squares solution of  Ax = b\n')
          print('The matrix A has %8g rows  and %8g cols' % (m, n))
          print('damp = %20.14e\n' % (damp))
          print('atol = %8.2e                 conlim = %8.2e\n' % (atol, conlim))
          print('btol = %8.2e             maxiter = %8g\n' % (btol, maxiter))
  
      u = b
      beta = norm(u)
  
      v = zeros(n)
      alpha = 0
  
      if beta > 0:
          u = (1 / beta) * u
          v = A.rmatvec(u)
          alpha = norm(v)
  
      if alpha > 0:
          v = (1 / alpha) * v
  
  
      # Initialize variables for 1st iteration.
  
      itn = 0
      zetabar = alpha * beta
      alphabar = alpha
      rho = 1
      rhobar = 1
      cbar = 1
      sbar = 0
  
      h = v.copy()
      hbar = zeros(n)
      x = zeros(n)
  
      # Initialize variables for estimation of ||r||.
  
      betadd = beta
      betad = 0
      rhodold = 1
      tautildeold = 0
      thetatilde = 0
      zeta = 0
      d = 0
  
      # Initialize variables for estimation of ||A|| and cond(A)
  
      normA2 = alpha * alpha
      maxrbar = 0
      minrbar = 1e+100
      normA = sqrt(normA2)
      condA = 1
      normx = 0
  
      # Items for use in stopping rules.
      normb = beta
      istop = 0
      ctol = 0
      if conlim > 0:
          ctol = 1 / conlim
      normr = beta
  
      # Reverse the order here from the original matlab code because
      # there was an error on return when arnorm==0
      normar = alpha * beta
      if normar == 0:
          if show:
              print(msg[0])
          return x, istop, itn, normr, normar, normA, condA, normx
  
      if show:
          print(' ')
          print(hdg1, hdg2)
          test1 = 1
          test2  = alpha / beta
          str1 = '%6g %12.5e' % (itn, x[0])
          str2 = ' %10.3e %10.3e' % (normr, normar)
          str3 = '  %8.1e %8.1e' % (test1,  test2)
          print(''.join([str1, str2, str3]))
  
      # Main iteration loop.
      while itn < maxiter:
          itn = itn + 1
  
          # Perform the next step of the bidiagonalization to obtain the
          # next  beta, u, alpha, v.  These satisfy the relations
          #         beta*u  =  a*v   -  alpha*u,
          #        alpha*v  =  A'*u  -  beta*v.
  
          u = A.matvec(v) - alpha * u
          beta = norm(u)
  
          if beta > 0:
              u = (1 / beta) * u
              v = A.rmatvec(u) - beta * v
              alpha = norm(v)
              if alpha > 0:
                  v = (1 / alpha) * v
  
          # At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
  
          # Construct rotation Qhat_{k,2k+1}.
  
          chat, shat, alphahat = _sym_ortho(alphabar, damp)
  
          # Use a plane rotation (Q_i) to turn B_i to R_i
  
          rhoold = rho
          c, s, rho = _sym_ortho(alphahat, beta)
          thetanew = s*alpha
          alphabar = c*alpha
  
          # Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar
  
          rhobarold = rhobar
          zetaold = zeta
          thetabar = sbar * rho
          rhotemp = cbar * rho
          cbar, sbar, rhobar = _sym_ortho(cbar * rho, thetanew)
          zeta = cbar * zetabar
          zetabar = - sbar * zetabar
  
          # Update h, h_hat, x.
  
          hbar = h - (thetabar * rho / (rhoold * rhobarold)) * hbar
          x = x + (zeta / (rho * rhobar)) * hbar
          h = v - (thetanew / rho) * h
  
          # Estimate of ||r||.
  
          # Apply rotation Qhat_{k,2k+1}.
          betaacute = chat * betadd
          betacheck = -shat * betadd
  
          # Apply rotation Q_{k,k+1}.
          betahat = c * betaacute
          betadd = -s * betaacute
  
          # Apply rotation Qtilde_{k-1}.
          # betad = betad_{k-1} here.
  
          thetatildeold = thetatilde
          ctildeold, stildeold, rhotildeold = _sym_ortho(rhodold, thetabar)
          thetatilde = stildeold* rhobar
          rhodold = ctildeold * rhobar
          betad = - stildeold * betad + ctildeold * betahat
  
          # betad   = betad_k here.
          # rhodold = rhod_k  here.
  
          tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold
          taud = (zeta - thetatilde * tautildeold) / rhodold
          d = d + betacheck * betacheck
          normr = sqrt(d + (betad - taud)**2 + betadd * betadd)
  
          # Estimate ||A||.
          normA2 = normA2 + beta * beta
          normA = sqrt(normA2)
          normA2 = normA2 + alpha * alpha
  
          # Estimate cond(A).
          maxrbar = max(maxrbar, rhobarold)
          if itn > 1:
            minrbar= min(minrbar, rhobarold)
          condA = max(maxrbar, rhotemp) / min(minrbar, rhotemp)
  
          # Test for convergence.
  
          # Compute norms for convergence testing.
          normar = abs(zetabar)
          normx = norm(x)
  
          # Now use these norms to estimate certain other quantities,
          # some of which will be small near a solution.
  
          test1 = normr / normb
          if (normA * normr) != 0:
              test2 = normar / (normA * normr)
          else:
              test2 = infty
          test3 = 1 / condA
          t1 = test1 / (1 + normA * normx / normb)
          rtol = btol + atol * normA * normx / normb
  
          # The following tests guard against extremely small values of
          # atol, btol or ctol.  (The user may have set any or all of
          # the parameters atol, btol, conlim  to 0.)
          # The effect is equivalent to the normAl tests using
          # atol = eps,  btol = eps,  conlim = 1/eps.
  
          if itn >= maxiter:
              istop = 7
          if 1 + test3 <= 1:
              istop = 6
          if 1 + test2 <= 1:
              istop = 5
          if 1 + t1 <= 1:
              istop = 4
  
          # Allow for tolerances set by the user.
  
          if test3 <= ctol:
              istop = 3
          if test2 <= atol:
              istop = 2
          if test1 <= rtol:
              istop = 1
  
          # See if it is time to print something.
  
          if show:
              if (n <= 40) or (itn <= 10) or (itn >= maxiter - 10) or \
                 (itn % 10 == 0) or (test3 <= 1.1 * ctol) or \
                 (test2 <= 1.1 * atol) or (test1 <= 1.1 * rtol) or \
                 (istop != 0):
  
                  if pcount >= pfreq:
                      pcount = 0
                      print(' ')
                      print(hdg1, hdg2)
                  pcount = pcount + 1
                  str1 = '%6g %12.5e' % (itn, x[0])
                  str2 = ' %10.3e %10.3e' % (normr, normar)
                  str3 = '  %8.1e %8.1e' % (test1, test2)
                  str4 = ' %8.1e %8.1e' % (normA, condA)
                  print(''.join([str1, str2, str3, str4]))
  
          if istop > 0:
              break
  
      # Print the stopping condition.
  
      if show:
          print(' ')
          print('LSMR finished')
          print(msg[istop])
          print('istop =%8g    normr =%8.1e' % (istop, normr))
          print('    normA =%8.1e    normAr =%8.1e' % (normA, normar))
          print('itn   =%8g    condA =%8.1e' % (itn, condA))
          print('    normx =%8.1e' % (normx))
          print(str1, str2)
          print(str3, str4)
  
      return x, istop, itn, normr, normar, normA, condA, normx
  
apis: scipy.sparse.linalg.interface.aslinearoperator numpy.linalg.norm numpy.zeros