title: Why are scipy's sparse solvers giving incorrect answers?
body: <p><strong>Note</strong>: I made a stupid mistake in the code I originally posted, as Warren Weckesser notes. When corrected, some of the solvers give correct answers, but others give NaN's or incorrect answers. I also forgot to include the runtime warnings in the output; they're there now. I've modified the question accordingly. I may be able to use the solvers that work, but I'd be much happier if I understood why the others are failing. </p>
  
  <p>I am trying to solve sparse systems of linear equations using one or more of the solvers found in scipy.sparse.linalg. In test cases, where the systems are small enough to solve directly, some of the sparse solvers are giving incorrect answers, as seen in the following example:</p>
  
  <pre><code>import numpy as np
  import scipy.sparse as ss
  
  
  A = np.matrix([[ 0.,  0.,  0.,  0.,  0.,  1., -1., -0., -0., -0., -0.],
         [ 0.,  0.,  0.,  0.,  0.,  2., -0., -1., -0., -0., -0.],
         [ 0.,  0.,  0.,  0.,  0.,  2., -0., -0., -1., -0., -0.],
         [ 0.,  0.,  0.,  0.,  0.,  2., -0., -0., -0., -1., -0.],
         [ 0.,  0.,  0.,  0.,  0.,  1., -0., -0., -0., -0., -1.],
         [ 1.,  2.,  2.,  2.,  1.,  0., -0., -0., -0., -0., -0.],
         [-1.,  0.,  0.,  0.,  0.,  0., -1., -0., -0., -0., -0.],
         [ 0., -1.,  0.,  0.,  0.,  0., -0., -1., -0., -0., -0.],
         [ 0.,  0., -1.,  0.,  0.,  0., -0., -0., -1., -0., -0.],
         [ 0.,  0.,  0., -1.,  0.,  0., -0., -0., -0., -1., -0.],
         [ 0.,  0.,  0.,  0., -1.,  0., -0., -0., -0., -0., -1.]])
  b = np.matrix([0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.]).T
  As = ss.coo_matrix(A)
  
  # The linear system Ax = b has a solution:
  x1 = np.linalg.solve(A,b)
  print("Solution to Ax = b:",x1)
  print("Ax - b = ",A*x1-b)
  
  print("Info and maximum error in solutions found by various other methods: ")
  x2,info = ss.linalg.bicg(As,b)
  print("bicg:",info,np.max(np.abs(x2-x1.ravel())))
  x2,info = ss.linalg.bicgstab(As,b)
  print("bicgstab:",info,np.max(np.abs(x2-x1.ravel())))
  x2,info = ss.linalg.cgs(As,b)
  print("cgs:",info,np.max(np.abs(x2-x1.ravel())))
  x2,info = ss.linalg.gmres(As,b)
  print("gmres:",info,np.max(np.abs(x2-x1.ravel())))
  x2,info = ss.linalg.lgmres(As,b)
  print("lgmres:",info,np.max(np.abs(x2-x1.ravel())))
  x2,info = ss.linalg.minres(As,b)
  print("minres:",info,np.max(np.abs(x2-x1.ravel())))
  x2,info = ss.linalg.qmr(As,b)
  print("qmr:",info,np.max(np.abs(x2-x1.ravel())))
  </code></pre>
  
  <p>When I run this, I get the following output:</p>
  
  <pre><code>Solution to Ax = b: [[ 0.07142857]
   [ 0.14285714]
   [ 0.14285714]
   [ 0.14285714]
   [ 0.07142857]
   [-0.07142857]
   [-0.07142857]
   [-0.14285714]
   [-0.14285714]
   [-0.14285714]
   [-0.07142857]]
  Ax - b =  [[  0.00000000e+00]
   [  0.00000000e+00]
   [  2.77555756e-17]
   [  0.00000000e+00]
   [  1.38777878e-17]
   [  0.00000000e+00]
   [  2.77555756e-17]
   [ -2.77555756e-17]
   [ -5.55111512e-17]
   [  2.77555756e-17]
   [  0.00000000e+00]]
  Info and maximum error in solutions found by various other methods: 
  bicg: 1 nan
  bicgstab: 1 nan
  cgs: 1 nan
  gmres: 0 5.55111512313e-17
  lgmres: 0 1.38777878078e-16
  minres: 0 0.142857142857
  qmr: -11 0.142857142857
  /Users/ebunn/anaconda/lib/python3.5/site-packages/scipy/sparse/linalg/isolve/iterative.py:197: RuntimeWarning: invalid value encountered in multiply
    work[slice2] *= sclr2
  /Users/ebunn/anaconda/lib/python3.5/site-packages/scipy/sparse/linalg/isolve/iterative.py:318: RuntimeWarning: invalid value encountered in multiply
    work[slice2] *= sclr2
  /Users/ebunn/anaconda/lib/python3.5/site-packages/scipy/sparse/linalg/isolve/minres.py:244: RuntimeWarning: divide by zero encountered in double_scalars
    Acond = gmax/gmin
  </code></pre>
  
  <p>Each calculation of x2 is supposed to be a solution to the same linear system as x1, so all of those errors in the last seven lines should be zero (or at least small).</p>
  
  <p>gmres and lgmres work, but the others don't. In most cases, the info correctly indicates failure, but minres indicates success (info=0) while returning a solution that is all zeros (incorrect).</p>
  
  <p>Here is some potentially relevant additional information:</p>
  
  <ul>
  <li>The matrix A is symmetric but not positive definite.</li>
  <li>The matrix A is quite well-conditioned.</li>
  <li>The results are unchanged if I replace the sparse representation As by the original matrix A in all of the sparse solvers -- that is, if I say <em>ss.linalg.bicg(A,b)</em> instead of <em>ss.linalg.bicg(As,b)</em>, etc.</li>
  <li>It's actually not fair to include bicgstab, because the documentation says that that is intended only for positive definite matrices. I included it in the hope that it might work, because this method is in principle applicable to indefinite matrices.</li>
  <li>I'm running Python 3.5.1, scipy version 0.17.0.</li>
  </ul>
  
  <p>Of course, for this system, it doesn't matter, because it's trivial to solve it directly. This is a warm-up problem for larger problems where sparsity will be essential. </p>
  
  <p>It could be that gmres and/or lgmres will meet my needs, but I'd still like to understand what's going wrong with the others, partly for my own understanding but also so that I might have a broader set of methods to choose from. In particular, neither of the two methods that work take advantage of A's symmetry, and it might be nice to have a method that does.</p>