Repository: scipy/scipy
Pull request: BUG: SuperLU 'NATURAL' order applies a column permutation
When using perm_c="NATURAL" in splu and spilu, the column ordering of a SuperLU
factorization is not preserved. This patch sets appropriate options to pass
into SuperLU so that the column ordering is preserved.

Row permutations for partial pivoting can still be performed.

Fix #7700, #8795

Problem statement:
SuperLU does not respect perm_c="NATURAL" 
According to 2.5.1 of the SuperLU user's guide, the `ColPerm` option (set via `scipy.sparse.linalg.splu`'s `permc_spec` parameter) should accept the value:
`NATURAL`: use natural ordering (i.e., Pc = I).
According to this description, I would expect that the returned SuperLU object's `perm_c` array to list the indices of the columns in increasing order (i.e. identical to `np.arange(m)`, where `m` is the number of columns). Sometimes this is the case, but it does not hard to find situations for which this option does not appear to be respected. Note that the order produced is different from the produced using the other options; but it is not the "natural" ordering.

Really, I'd like to use `permc_spec="MY_PERMC"` , which appears to be a valid value, but the parameter for actually specifying `perm_c` does not appear to be exposed.
 
### Reproducing code example:
```python
np.random.seed(1)
m = 1000
p = 0.01
A = scipy.sparse.random(m,m,p)
f = splu(A, permc_spec = "NATURAL")
# f = splu(A, options={"ColPerm":"NATURAL"}) # no difference
print(np.allclose(f.perm_c,np.arange(m)))
print(f.perm_c[:10])
```

### Error message:
Output from code above:
```
False
[21 22  5 17  8 23 15 24  6 20]
```

### Scipy/Numpy/Python version information:
```
('0.19.0', '1.11.3', sys.version_info(major=2, minor=7, micro=11, releaselevel='final', serial=0))
```

Hints from issue discussion:
I have been working with SuperLU (and sparse matrices in general) quite a bit lately and I have many suggestions, one of which would be to remove the restriction:
```python
M, N = A.shape
if (M != N):
    raise ValueError("can only factor square matrices")  # is this true?
```
as the first paragraph of the introduction to the SuperLU user's guide states:
"The LU factorization routines can handle non-square matrices."
Should I propose these as enhancements?
Should I refrain from suggestions about sparse matrices themselves? (like add the sparse equivalent for `np.allclose()`? I get the impression that there are enough known issues.) 


Suggestions for enhancements are welcome (as always, with the understanding that it all runs on volunteer effort). Please file separate tickets, though.
The arguments are parsed correctly. The question is then reading this code here:
https://github.com/scipy/scipy/blob/ff052f6e7c5f952f93ebceeec58159e64693e806/scipy/sparse/linalg/dsolve/_superluobject.c#L737
and comparing with the SuperLU manual. It's not clear it does the expected thing.

The input parameters with just the ColPerm option are
```
options  = {Fact = DOFACT, Equil = YES, ColPerm = NATURAL, Trans = NOTRANS, IterRefine = NOREFINE, 
  DiagPivotThresh = 1, SymmetricMode = NO, PivotGrowth = NO, ConditionNumber = NO, 
  RowPerm = NOROWPERM, ILU_DropRule = 1, ILU_DropTol = 0, 
  ILU_FillFactor = 4.6355710338275721e-310, ILU_Norm = (unknown: 10000), 
  ILU_FillTol = 6.9533558071986478e-310, ILU_MILU = 1000, 
  ILU_MILU_Dim = 6.9533558071946952e-310, ParSymbFact = (YES | unknown: 3859408010), 
  ReplaceTinyPivot = (YES | unknown: 32766), SolveInitialized = (unknown: 1441640144), 
  RefineInitialized = (YES | unknown: 21844), PrintStat = YES, nnzL = 21845, nnzU = 0, 
  num_lookaheads = 0, lookahead_etree = YES, SymPattern = NO}
```
I have a fix to this (workaround, really). Simply call
```python
f = splu(A, permc_spec = "NATURAL", diag_pivot_thresh=0, options={"SymmetricMode":True})
```
From 2.6 of SuperLU docs (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ug.pdf), `SymmetricMode` "prefers" the diagonal pivot, but since it's not specific as to how strongly it prefers this, setting `diag_pivot_thresh=0` as well should guarantee that no pivoting is done.

(Note it is strange that setting `diag_pivot_thresh=0` alone is not enough)

Commit message: TST: linalg: improve SuperLU permc_spec=NATURAL tests

Patch to review:
diff --git a/scipy/sparse/linalg/dsolve/linsolve.py b/scipy/sparse/linalg/dsolve/linsolve.py
index f12ba71fc7b8..93f908c00b6f 100644
--- a/scipy/sparse/linalg/dsolve/linsolve.py
+++ b/scipy/sparse/linalg/dsolve/linsolve.py
@@ -321,6 +321,11 @@ def splu(A, permc_spec=None, diag_pivot_thresh=None,
                     PanelSize=panel_size, Relax=relax)
     if options is not None:
         _options.update(options)
+
+    # Ensure no column permutations applied
+    if (_options["ColPerm"] == "NATURAL"):
+        _options["SymmetricMode"] = True
+
     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
                           csc_construct_func=csc_construct_func,
                           ilu=False, options=_options)
@@ -407,6 +412,11 @@ def spilu(A, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None,
                     PanelSize=panel_size, Relax=relax)
     if options is not None:
         _options.update(options)
+
+    # Ensure no column permutations applied
+    if (_options["ColPerm"] == "NATURAL"):
+        _options["SymmetricMode"] = True
+
     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
                           csc_construct_func=csc_construct_func,
                           ilu=True, options=_options)
diff --git a/scipy/sparse/linalg/dsolve/tests/test_linsolve.py b/scipy/sparse/linalg/dsolve/tests/test_linsolve.py
index f8cb99092da0..ada45ffc8fbb 100644
--- a/scipy/sparse/linalg/dsolve/tests/test_linsolve.py
+++ b/scipy/sparse/linalg/dsolve/tests/test_linsolve.py
@@ -547,6 +547,31 @@ def test_splu_perm(self):
         lu = splu(a_)
         assert_array_equal(lu.perm_r, lu.perm_c)
 
+    @pytest.mark.parametrize("splu_fun, rtol", [(splu, 1e-7), (spilu, 1e-1)])
+    def test_natural_permc(self, splu_fun, rtol):
+        # Test that the "NATURAL" permc_spec does not permute the matrix
+        np.random.seed(42)
+        n = 500
+        p = 0.01
+        A = scipy.sparse.random(n, n, p)
+        x = np.random.rand(n)
+        # Make A diagonal dominant to make sure it is not singular
+        A += (n+1)*scipy.sparse.identity(n)
+        A_ = csc_matrix(A)
+        b = A_ @ x
+
+        # without permc_spec, permutation is not identity
+        lu = splu_fun(A_)
+        assert_(np.any(lu.perm_c != np.arange(n)))
+
+        # with permc_spec="NATURAL", permutation is identity
+        lu = splu_fun(A_, permc_spec="NATURAL")
+        assert_array_equal(lu.perm_c, np.arange(n))
+
+        # Also, lu decomposition is valid
+        x2 = lu.solve(b)
+        assert_allclose(x, x2, rtol=rtol)
+
     @pytest.mark.skipif(not hasattr(sys, 'getrefcount'), reason="no sys.getrefcount")
     def test_lu_refcount(self):
         # Test that we are keeping track of the reference count with splu.


Code review comments:
diff_hunk: @@ -407,6 +412,11 @@ def spilu(A, drop_tol=None, fill_factor=None, drop_rule=None, permc_spec=None,
                     PanelSize=panel_size, Relax=relax)
     if options is not None:
         _options.update(options)
+
+    # Ensure no column permutations applied
path: scipy/sparse/linalg/dsolve/linsolve.py
text: @user1:
```suggestion
    # Ensure that no column permutations are applied
```

diff_hunk: @@ -321,6 +321,11 @@ def splu(A, permc_spec=None, diag_pivot_thresh=None,
                     PanelSize=panel_size, Relax=relax)
     if options is not None:
         _options.update(options)
+
+    # Ensure no column permutations applied
path: scipy/sparse/linalg/dsolve/linsolve.py
text: @user1:
```suggestion
    # Ensure that no column permutations are applied
```