BUG: sparse/dsolve: fix dense matrix fortran order bugs in superlu wrappers #3367

Merged
merged 4 commits into from Feb 23, 2014

3 participants

@pv
SciPy member
pv commented Feb 21, 2014

The splu() SuperLU wrapper has been broken since time immemorial (< 0.7.2)
for matrix inputs, because of confusion in the C code between Fortran and C
order inputs.

This is fairly simple to fix, so let's fix it. This also allows passing in dense matrices
directly to SuperLU, without needing to loop on the Python side.

Fixes gh-3363

@pv pv BUG: sparse/dsolve: fix dense matrix fortran order bugs in superlu wr…
…appers

This also allows passing in dense matrices to SuperLU, without needing
to loop on the Python side.
2344e41
@pv pv added this to the 0.14.0 milestone Feb 21, 2014
@pv pv added PR and removed prio-high labels Feb 21, 2014
@coveralls

Coverage Status

Changes Unknown when pulling 2344e41 on pv:splu-bugfix into ** on scipy:master**.

@coveralls

Coverage Status

Coverage remained the same when pulling a18a98a on pv:splu-bugfix into 1b34692 on scipy:master.

@rgommers rgommers commented on the diff Feb 23, 2014
scipy/sparse/linalg/dsolve/_superluobject.c
-solves linear system of equations with one or sereral right hand sides.\n\
-\n\
-parameters\n\
-----------\n\
-\n\
-b array, right hand side(s) of equation\n\
-x array, solution vector(s)\n\
-trans 'N': solve A * x == b\n\
- 'T': solve A^T * x == b\n\
- 'H': solve A^H * x == b\n\
- (optional, default value 'N')\n\
-";
+static char solve_doc[] = (
+ "x = self.solve(b, trans)\n"
+ "\n"
+ "solves linear system of equations with one or sereral right hand sides.\n"
@rgommers
SciPy member
rgommers added a line comment Feb 23, 2014

Might as well fix the typo here (sereral).

@pv
SciPy member
pv added a line comment Feb 23, 2014

I fix it in gh-3375 :)

@rgommers
SciPy member
rgommers added a line comment Feb 23, 2014

OK

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers and 1 other commented on an outdated diff Feb 23, 2014
scipy/sparse/linalg/dsolve/linsolve.py
else:
- if isspmatrix_csc(A):
- flag = 1 # CSC format
- else:
- flag = 0 # CSR format
- options = dict(ColPerm=permc_spec)
- x = _superlu.gssv(N, A.nnz, A.data, A.indices, A.indptr, b_vec, flag,
- options=options)[0]
- else:
- # Cover the case where b is a matrix
- if b_is_sparse:
+ flag = 0 # CSR format
@rgommers
SciPy member
rgommers added a line comment Feb 23, 2014

flag isn't used in the else clause below, so why distinguish here between CSC and CSR? Something is missing below, or this check can be removed.

@pv
SciPy member
pv added a line comment Feb 23, 2014

flag is used on line 144 below

@rgommers
SciPy member
rgommers added a line comment Feb 23, 2014

Yes, but that's in the code handling non-sparse b, so shouldn't care about CSR vs CSC. And if it's only used there, can be moved in the if branch.

@pv
SciPy member
pv added a line comment Feb 23, 2014

Ok, moved it down inside the if clause.

@rgommers
SciPy member
rgommers added a line comment Feb 23, 2014

Sorry if I'm missing something obvious, but if b_is_sparse is False, then checking with isspmatrix_csc is still useless no? flag can only be 0.

@pv
SciPy member
pv added a line comment Feb 23, 2014

It's checking the matrix format of A, not that of b.
A can be either CSC or CSR.

@rgommers
SciPy member
rgommers added a line comment Feb 23, 2014

OK I'm being dense.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
SciPy member

Doesn't build with Bento (edit: also not with distutils):

[ 289/1158] c: scipy/optimize/Zeros/brenth.c -> build/scipy/optimize/Zeros/brenth.c.10.o
In file included from /home/rgommers/.local/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:14:0,
                 from ../scipy/sparse/linalg/dsolve/_superluobject.h:19,
                 from ../scipy/sparse/linalg/dsolve/_superluobject.c:13:
../scipy/sparse/linalg/dsolve/_superluobject.c: In function ‘SciPyLU_solve’:
../scipy/sparse/linalg/dsolve/_superluobject.c:85:9: error: ‘NPY_ARRAY_F_CONTIGUOUS’ undeclared (first use in this function)
         NPY_ARRAY_F_CONTIGUOUS | NPY_ARRAY_ENSURECOPY);
         ^
/home/rgommers/.local/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:98:28: note: in definition of macro ‘PyArray_FROMANY’
                         (((flags) & NPY_ENSURECOPY) ?                         \
                            ^
compilation terminated due to -Wfatal-errors.
In file included from ../scipy/sparse/sparsetools/py3k.h:23:0,
                 from ../scipy/sparse/sparsetools/coo_wrap.cxx:3054:
/home/rgommers/.local/lib/python2.7/site-packages/numpy/core/include/numpy/npy_3kcompat.h: In function ‘PyObject* npy_PyFile_OpenFile(PyObject*, char*)’:
/home/rgommers/.local/lib/python2.7/site-packages/numpy/core/include/numpy/npy_3kcompat.h:258:60: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
     return PyObject_CallFunction(open, "Os", filename, mode);
@rgommers
SciPy member

Since this builds on Travis, likely related to numpy version (I use 1.5.1).

@pv
SciPy member
pv commented Feb 23, 2014

The failure is due to Numpy 1.5. Fixed. (We should make one of the travis builds use this Numpy version, so these issues would be easier to catch.)

@coveralls

Coverage Status

Coverage remained the same when pulling 8b9cdef on pv:splu-bugfix into 1b34692 on scipy:master.

@rgommers
SciPy member

Works for me now, so in it goes. Thanks @pv

@rgommers rgommers merged commit 25125e5 into scipy:master Feb 23, 2014

1 check passed

Details default The Travis CI build passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment