Gmres with Sparse Complex Matrices gives Segmentation fault (or NULL result)

Hi all, I am having some troubles using iterative solvers with sparse complex matrices. (I first posted that issue on the SciPy-User mailing list but I guess this one is more appropriate.) So, here is an example with gmres: the script from scipy import sparse from numpy.linalg import norm from numpy.random import rand import scipy.sparse.linalg as spla C = sparse.lil_matrix((10, 10), dtype=complex) C.setdiag(rand(10)) C[0,3] = 1j C = C.tocsr() c = rand(10)+rand(10)*1j x = spla.spsolve(C, c) print "spsolve: norm(C*x - c) = ", norm(C*x - c) (x,info) = spla.gmres(C, c) print "gmres: norm(C*x - c) = ", norm(C*x - c) gives as output: spsolve: norm(C*x - c) = 1.57621370006e-16 Segmentation fault Actually, sometimes I get this error message instead of a Segmentation fault: Traceback (most recent call last): File "test_gmres_cplx.py", line 45, in <module> (x,info) = spla.gmres(C, c) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/iterative.py", line 437, in gmres revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob) SystemError: NULL result without error in PyObject_Call Note that I have tested gmres with real sparse matrices, and it runs fine: Indeed C = sparse.lil_matrix((10, 10)) C.setdiag(rand(10)) C = C.tocsr() c = rand(10) x = spla.spsolve(C, c) print "spsolve: norm(C*x - c) = ", norm(C*x - c) (x,info) = spla.gmres(C, c) print "gmres: norm(C*x - c) = ", norm(C*x - c) gives spsolve: norm(C*x - c) = 6.93889390391e-18 gmres: norm(C*x - c) = 5.28860261481e-16 I have looked on the web for solutions but haven't found any. Some very old posts indicate similar errors but they don't come with an answer, and I imagine that if those were due to bugs, they should have been fixed by now... Am I doing something stupid here, or is that a real problem ? Is somebody aware of a solution ? (I am using scipy version 0.10.1) Thanks, Martin -- Martin Campos Pinto LJLL, UPMC & CNRS

24.11.2012 07:50, Martin Campos Pinto kirjoitti: [clip]
I have looked on the web for solutions but haven't found any. Some very old posts indicate similar errors but they don't come with an answer, and I imagine that if those were due to bugs, they should have been fixed by now...
Am I doing something stupid here, or is that a real problem ? Is somebody aware of a solution ? (I am using scipy version 0.10.1)
I don't get a crash with this code example. Valgrind also does not show anything, so I suppose this is something occurring only on OSX. There are some known issues on OSX with Accelerate, but I'm not sure if they are in play here. Can you run `scipy.test("full", verbose=2)` without seeing problems? -- Pauli Virtanen

Hi, thank you for the quick answer. The run `scipy.test("full", verbose=2)' fails with the following error message (I'm only pasting the last few lines, but I can send the whole output if necessary...) ====================================================================== ERROR: Failure: ImportError (dlopen(/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/atlas_version.so, 2): Symbol not found: _ATL_buildinfo Referenced from: /Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/atlas_version.so Expected in: dynamic lookup ) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/nose/loader.py", line 390, in loadTestsFromName addr.filename, addr.module) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/nose/importer.py", line 39, in importFromPath return self.importFromDir(dir_path, fqname) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/nose/importer.py", line 86, in importFromDir mod = load_module(part_fqname, fh, filename, desc) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/tests/test_atlas_version.py", line 6, in <module> import scipy.linalg.atlas_version ImportError: dlopen(/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/atlas_version.so, 2): Symbol not found: _ATL_buildinfo Referenced from: /Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/atlas_version.so Expected in: dynamic lookup ====================================================================== FAIL: test_dot (test_blas.TestFBLAS1Simple) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/lib/blas/tests/test_blas.py", line 71, in test_dot assert_almost_equal(f([3j,-4,3-4j],[2,3,1]),-9+2j) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/numpy/testing/utils.py", line 448, in assert_almost_equal raise AssertionError(msg) AssertionError: Arrays are not almost equal to 7 decimals ACTUAL: 0j DESIRED: (-9+2j) ====================================================================== FAIL: test_complex_dotc (test_blas.TestFBLAS1Simple) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/tests/test_blas.py", line 121, in test_complex_dotc assert_almost_equal(f([3j,-4,3-4j],[2,3j,1]),3-14j) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/numpy/testing/utils.py", line 448, in assert_almost_equal raise AssertionError(msg) AssertionError: Arrays are not almost equal to 7 decimals ACTUAL: 0j DESIRED: (3-14j) ====================================================================== FAIL: test_complex_dotu (test_blas.TestFBLAS1Simple) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/linalg/tests/test_blas.py", line 115, in test_complex_dotu assert_almost_equal(f([3j,-4,3-4j],[2,3,1]),-9+2j) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/numpy/testing/utils.py", line 448, in assert_almost_equal raise AssertionError(msg) AssertionError: Arrays are not almost equal to 7 decimals ACTUAL: 0j DESIRED: (-9+2j) ---------------------------------------------------------------------- Ran 5045 tests in 899.955s FAILED (KNOWNFAIL=14, SKIP=29, errors=1, failures=3) Le 11/24/12 10:40 AM, Pauli Virtanen a écrit :
24.11.2012 07:50, Martin Campos Pinto kirjoitti: [clip]
I have looked on the web for solutions but haven't found any. Some very old posts indicate similar errors but they don't come with an answer, and I imagine that if those were due to bugs, they should have been fixed by now...
Am I doing something stupid here, or is that a real problem ? Is somebody aware of a solution ? (I am using scipy version 0.10.1) I don't get a crash with this code example. Valgrind also does not show anything, so I suppose this is something occurring only on OSX.
There are some known issues on OSX with Accelerate, but I'm not sure if they are in play here.
Can you run `scipy.test("full", verbose=2)` without seeing problems?
-- Martin Campos Pinto http://www-irma.u-strasbg.fr/~campos

Hi, I have just tried with scipy version 0.11.0, installed with mac ports and gmres runs without crashing. Since scipy 0.10.1 was the package provided with EPD I first thought it should run fine... Thank you, Martin -- ps. Also, just to check I have run your test `scipy.test("full", verbose=2)' with scipy 0.11.0 and there still is a problem. Python 2.7.3 (default, Oct 22 2012, 06:12:28) [GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin Type "help", "copyright", "credits" or "license" for more information.
import scipy scipy.test("full", verbose=2) Running unit tests for scipy NumPy version 1.6.2 NumPy is installed in /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy SciPy version 0.11.0 SciPy is installed in /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy Python version 2.7.3 (default, Oct 22 2012, 06:12:28) [GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] nose version 1.2.1 Tests cophenet(Z) on tdist data set. ... ok Tests cophenet(Z, Y) on tdist data set. ... ok Tests correspond(Z, y) on linkage and CDMs over observation sets of different sizes. ... ok Tests correspond(Z, y) on linkage and CDMs over observation sets of different sizes. Correspondance should be false. ... ok Tests correspond(Z, y) on linkage and CDMs over observation sets of different sizes. Correspondance should be false. ... ok Tests correspond(Z, y) with empty linkage and condensed distance matrix. ... ok Tests num_obs_linkage with observation matrices of multiple sizes. ... ok Tests dendrogram calculation on single linkage of the tdist data set. ... ok Tests fcluster(Z, criterion='maxclust', t=2) on a random 3-cluster data set. ... ok
(...) test_exclusive_end (test_slice_handler.TestBuildSliceAtom) ... ok match slice from a[1:] ... ok match slice from a[1::] ... ok match slice from a[1:2] ... ok match slice from a[1:2:] ... ok match slice from a[1:2:3] ... ok match slice from a[1::3] ... ok match slice from a[:] ... ok match slice from a[::] ... ok match slice from a[:2] ... ok match slice from a[:2:] ... ok match slice from a[:2:3] ... ok match slice from a[:1+i+2:] ... ok match slice from a[0] ... ok match slice from a[::3] ... ok transform a[:] to slice notation ... ok transform a[:,:] = b[:,1:1+2:3] *(c[1-2+i:,:] - c[:,:]) ... ok test_type_match_array (test_standard_array_spec.TestArrayConverter) ... ok test_type_match_int (test_standard_array_spec.TestArrayConverter) ... ok test_type_match_string (test_standard_array_spec.TestArrayConverter) ... ok ====================================================================== FAIL: test_iv_cephes_vs_amos_mass_test (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/special/tests/test_basic.py", line 1659, in test_iv_cephes_vs_amos_mass_test assert_(dc[k] < 2e-7, (v[k], x[k], special.iv(v[k], x[k]), special.iv(v[k], x[k]+0j))) File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/testing/utils.py", line 34, in assert_ raise AssertionError(msg) AssertionError: (189.2947429454936, 3.0238805556481037, 4.089165443940765e-317, 0j) ---------------------------------------------------------------------- Ran 6217 tests in 1046.384s FAILED (KNOWNFAIL=15, SKIP=32, failures=1) <nose.result.TextTestResult run=6217 errors=0 failures=1> Le 11/24/12 10:40 AM, Pauli Virtanen a écrit :
24.11.2012 07:50, Martin Campos Pinto kirjoitti: [clip]
I have looked on the web for solutions but haven't found any. Some very old posts indicate similar errors but they don't come with an answer, and I imagine that if those were due to bugs, they should have been fixed by now...
Am I doing something stupid here, or is that a real problem ? Is somebody aware of a solution ? (I am using scipy version 0.10.1) I don't get a crash with this code example. Valgrind also does not show anything, so I suppose this is something occurring only on OSX.
There are some known issues on OSX with Accelerate, but I'm not sure if they are in play here.
Can you run `scipy.test("full", verbose=2)` without seeing problems?
-- Martin Campos Pinto http://www-irma.u-strasbg.fr/~campos

Hi, 28.11.2012 06:56, Martin Campos Pinto kirjoitti:
I have just tried with scipy version 0.11.0, installed with mac ports and gmres runs without crashing. Since scipy 0.10.1 was the package provided with EPD I first thought it should run fine...
That is sort of strange in itself, given that there are essentially no changes between 0.10.1 and 0.11.0 in that part of code. The test suite also includes complex-valued test problems for gmres, so it not crashing is also not immediately understood. If the crashing version was 0.9.0 or earlier, then I'd understand what is going on... If you still can, double-checking what 'print scipy.__version__' said could be useful. But nice to hear it works now. (And for future reference, PyAMG contains pure-python Krylov solvers, including gmres.) Best, Pauli

Hi, Le 11/28/12 7:26 PM, Pauli Virtanen a écrit :
Hi,
28.11.2012 06:56, Martin Campos Pinto kirjoitti:
I have just tried with scipy version 0.11.0, installed with mac ports and gmres runs without crashing. Since scipy 0.10.1 was the package provided with EPD I first thought it should run fine... That is sort of strange in itself, given that there are essentially no changes between 0.10.1 and 0.11.0 in that part of code. The test suite also includes complex-valued test problems for gmres, so it not crashing is also not immediately understood.
If the crashing version was 0.9.0 or earlier, then I'd understand what is going on... If you still can, double-checking what 'print scipy.__version__' said could be useful.
I have just reinstalled my os (for reasons not related to scipy...) and updated the scipy mac ports install but not the EPD one, but that one already used scipy 0.10.0 and the result is the same. Here's the test with the scipy version number written down: import scipy from scipy import sparse from numpy.linalg import norm from numpy.random import rand import scipy.sparse.linalg as spla print "scipy version: ", scipy.__version__ C = sparse.lil_matrix((10, 10), dtype=complex) C.setdiag(rand(10)) C[0,3] = 1j C = C.tocsr() c = rand(10)+rand(10)*1j x = spla.spsolve(C, c) print "spsolve: norm(C*x - c) = ", norm(C*x - c) (x,info) = spla.gmres(C, c) print "gmres: norm(C*x - c) = ", norm(C*x - c) gives: 1. with the python binary installed with mac ports: Python 2.7.3 (default, Oct 22 2012, 06:12:28) Martins-MacBook-Pro:~ campos$ /opt/local/bin/python2.7 test_gmres_cplx.py scipy version: 0.11.0 spsolve: norm(C*x - c) = 2.60370378581e-16 gmres: norm(C*x - c) = 1.21827900786e-06 2. with the python binary installed with EPD: Python 2.7.2 |EPD 7.2-2 (64-bit)| (default, Sep 7 2011, 16:31:15) Martins-MacBook-Pro:~ campos$ python test_gmres_cplx.py scipy version: 0.10.0 spsolve: norm(C*x - c) = 2.48641008715e-16 Segmentation fault Best, Martin -- Martin Campos Pinto http://www-irma.u-strasbg.fr/~campos
participants (2)
-
Martin Campos Pinto
-
Pauli Virtanen