Interest in adopting scipy patches?

[dschmitz89]

SuperLU is the workhorse for sparse linear solves in scipy. Currently, we vendor a copy of the superLU source and apply a few patches. Ideally, we would be able to just use your code as a submodule (or subproject so that a system installed version of superLU can be used) so that we can update without recreating the patch file etc (see discussion at scipy) . Therefore I would like to ask if our patches are of interest to you. The full patch file can be found here. SciPy is currently at commit 18958c3.

From what I understand, our patches are mostly about avoiding crashes in BLAS function for singular matrices. An example is the following:

--- a/scipy/sparse/linalg/_dsolve/SuperLU/SRC/cpanel_bmod.c
+++ b/scipy/sparse/linalg/_dsolve/SuperLU/SRC/cpanel_bmod.c
@@ -428,6 +428,12 @@ cpanel_bmod (
 		    CTRSV( ftcs1, ftcs2, ftcs3, &segsze, &lusup[luptr], 
 			   &nsupr, tempv, &incx );
 #else
+#if SCIPY_FIX
+		   if (nsupr < segsze) {
+			/* Fail early rather than passing in invalid parameters to TRSV. */
+			ABORT("failed to factorize matrix");
+		   }
+#endif

Would you accept a PR that applies such patches to SuperLU directly?

[dschmitz89]

CC @xiaoyeli

[xiaoyeli]

@dschmitz89 Sorry for the delayed response. I was traveling a lot.

The singluarity should be already handled in pivoting routine dpivotL.c, see the code here:

superlu/SRC/dpivotL.c

Line 135 in 87a1d1a

/* Test for singularity */

I suggest that in the caller dgstrf.c, after the call to dpivotL(), check if return value is not zero, then you abort the program, see code here:

superlu/SRC/dgstrf.c

Line 333 in 87a1d1a

if ( (*info = dpivotL(icol, diag_pivot_thresh, &usepr, perm_r,

superlu/SRC/dgstrf.c

Line 388 in 87a1d1a

if ( (*info = dpivotL(jj, diag_pivot_thresh, &usepr, perm_r,

This avoids going to all the downstream BLAS level to check.

If you set up a PR in this cleaner manner, I am happy to merge it. Thanks.

[xiaoyeli]

@dschmitz89 Checking to see whether what I proposed 3 weeks ago makes sense to you.

[dschmitz89]

@xiaoyeli Your proposal makes a lot of sense to me. I will see if I can find the time for a PR.

[gruenich]

@dschmitz89 Are you still interested in creating a PR to get rid of your scipy patches? It would be a shame to let this idea bit-rot.

[dschmitz89]

I didn't have the time so far. If anyone beats me to it, that's of course also fine.

[rgommers]

The singluarity should be already handled in pivoting routine dpivotL.c, see the code here:

I investigated this approach, and it doesn't work. The main reason is that the crashes come from this branch:

superlu/SRC/dgstrf.c

Line 345 in 53afbff

} else { /* Work on one panel of panel_size columns */

and the linked dpivotL call on line 333 is in the if branch, so it doesn't get hit. The second dpivotL call in that else branch only happens after dpanel_bmod is already called.

There are more LAPACK calls with invalid parameters with similar issues. For example, in the if branch the dsnode_bmod call has a problem, and it's also above the dpivotL call.

superlu/SRC/dgstrf.c

Line 331 in 53afbff

dsnode_bmod(icol, jsupno, fsupc, dense, tempv, Glu, stat);

I believe the SciPy fixes are about as clean as it can be done. They also pass the SuperLU test suite. I'll submit a PR with one commit per problem, so it's easier to see what the actual changes are.

The SuperLU test suite doesn't cover the issues. The SciPy one does; if I disable the SciPy fixes I see:

 ** On entry to DTRSV  parameter number  6 had an illegal value
 ** On entry to DGEMV  parameter number  2 had an illegal value
 ** On entry to DTRSV  parameter number  6 had an illegal value
 ** On entry to DGEMV  parameter number  2 had an illegal value

Those messages come from XERBLA; the actual behavior depends on the LAPACK implementation - Netlib will abort, OpenBLAS/MKL mostly return to let the user handle it, but the output may be in an undefined state - so this should always be avoided.

I didn't try to convert the SciPy tests covering these problems to the SuperLU test suite. I can if desired/needed to get our fixes upstreamed. It looks a bit fiddly, but doable. Here is an example:

https://github.com/scipy/scipy/blob/ee66af543ecb5af36c912455d56a6ed167a0eaf6/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py#L210-L227