Repository: scipy/scipy
Pull request: ENH: sparse.linalg: speed up `spsolve_triangular`
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes #16131.


#### What does this implement/fix?
<!--Please explain your changes.-->
#16131 reports that in `scipy.sparse.linalg`, the [`spsolve_triangular`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve_triangular.html) function is much slower than the [`spsolve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html) function when both are applied to a triangular system. This is unexpected because `spsolve_triangular` is the more specialized function. 

I found that `spsolve_triangular` is slower because it is entirely implemented in Python, while `spsolve` calls a compiled function from SuperLU. This PR suggests to use the SuperLU [`gstrs` (general sparse triangular solve)](https://github.com/scipy/scipy/blob/80fccfa2e1776170220cc07e1b2871325dd64ea5/scipy/sparse/linalg/_dsolve/SuperLU/SRC/dgstrs.c#L100) family of functions for `spsolve_triangular`. 

On my system, the [reproducing code example](https://github.com/scipy/scipy/issues/16131#issue-1228053542) prints e.g.
```
spsolve_triangular: 14.718287706375122 seconds
spsolve: 0.5277576446533203 seconds
```
With the changes suggested here, the example prints e.g.
```
<some-path>/scipy/build-install/lib/python3/dist-packages/scipy/sparse/linalg/_dsolve/linsolve.py:637: SparseEfficiencyWarning: CSC matrix format is required. Converting to CSC matrix.
  warn('CSC matrix format is required. Converting to CSC matrix.',
spsolve_triangular: 0.08869504928588867 seconds
spsolve: 0.36919116973876953 seconds
```

#### Additional information
<!--Any additional information you think is important.-->
This is my first PR in SciPy and I have never worked with SuperLU before. 


Problem statement:
BUG: spsolve_triangular is way slower than spsolve
### Describe your issue.

Hi! 
I used `spsolve_trianguar` in a project and noticed that is is really slow. Then I compared it to the generic `spsolve` method and this one is indeed a lot faster. I created a mini program that simply solves `Ax=b` for an identity matrix `A` and a zero vector `b`

For `n=1000000` the measured time (on a single CPU core) is:
```
spsolve_triangular: 21.094783067703247 seconds
spsolve: 0.34332990646362305 seconds
```
On my more advanced triangular matrix the time difference is also in this range. Note that it includes the decomposition steps done in `spsolve`.

### Reproducing Code Example

```python
from scipy.sparse.linalg import spsolve
from scipy.sparse.linalg import spsolve_triangular

import numpy as np

import scipy
import time

n = 1000000

# create a sparse identity matrix and zero vector
I = scipy.sparse.identity(n).tocsr()
b = np.zeros(n)

#compare spsolve_triangular with spsolve

start = time.time()
x = spsolve_triangular(I, b)
stop = time.time()

print( "spsolve_triangular:",stop-start,"seconds")

start = time.time()
x = spsolve(I, b)
stop = time.time()

print( "spsolve:",stop-start,"seconds")
```


### Error message

```shell
None. But the generic algorithm is much faster than the specialized one. This should not happen.
```


### SciPy/NumPy/Python version information

1.7.3 1.21.5 sys.version_info(major=3, minor=10, micro=4, releaselevel='final', serial=0)

Hints from issue discussion:
~~I think this is a duplicate of gh-10722. Do you agree @pjaap?~~
Never mind! This is the sparse version of that issue.

Commit message: ENH: spsolve_triangular raises LinAlgError

instead of warning when the matrix is singular.

Co-authored-by: Jake Bowhay <60778417+j-bowhay@users.noreply.github.com>

Patch to review:
diff --git a/scipy/sparse/linalg/_dsolve/_superlumodule.c b/scipy/sparse/linalg/_dsolve/_superlumodule.c
index 6480afdc1f24..31b74451366e 100644
--- a/scipy/sparse/linalg/_dsolve/_superlumodule.c
+++ b/scipy/sparse/linalg/_dsolve/_superlumodule.c
@@ -17,6 +17,7 @@
 #include <numpy/ndarrayobject.h>
 
 #include "_superluobject.h"
+#include "SuperLU/SRC/superlu_enum_consts.h"
 
 
 /*
@@ -263,6 +264,170 @@ static PyObject *Py_gstrf(PyObject * self, PyObject * args,
     return NULL;
 }
 
+static PyObject *Py_gstrs(PyObject * self, PyObject * args,
+                          PyObject * keywds)
+{
+  /* compressed sparse row matrix L */
+  int L_N=0, L_nnz=0;
+  PyArrayObject *L_nzvals=NULL, *L_colind=NULL, *L_rowptr=NULL;
+  /* compressed sparse column matrix U */
+  int U_N=0, U_nnz=0;
+  PyArrayObject *U_nzvals=NULL, *U_rowind=NULL, *U_colptr=NULL;
+  /* right hand side / solution */
+  PyObject *X_py=NULL;
+  /* whether the matrix is transposed ('T'), conjugate transposed ('H') or normal ('N') */
+  volatile int itrans = 'N';
+  volatile jmp_buf* jmpbuf_ptr;
+  SLU_BEGIN_THREADS_DEF;
+
+
+  static char* kwlist[] = {"trans",
+                           "L_N", "L_nnz", "L_nzvals", "L_colind", "L_rowptr",
+                           "U_N", "U_nnz", "U_nzvals", "U_rowind", "U_colptr",
+                           "B", NULL};
+
+  /* Parse and check input arguments. */
+
+  int res =
+      PyArg_ParseTupleAndKeywords(args, keywds, "CiiO!O!O!iiO!O!O!O", kwlist,
+        &itrans,
+        &L_N, &L_nnz, &PyArray_Type, &L_nzvals, &PyArray_Type, &L_colind, &PyArray_Type, &L_rowptr,
+        &U_N, &U_nnz, &PyArray_Type, &U_nzvals, &PyArray_Type, &U_rowind, &PyArray_Type, &U_colptr,
+        &X_py );
+  if(!res)
+    return NULL;
+
+  volatile trans_t trans; 
+  if (itrans == 'n' || itrans == 'N')
+      trans = NOTRANS;
+  else if (itrans == 't' || itrans == 'T')
+      trans = TRANS;
+  else if (itrans == 'h' || itrans == 'H')
+      trans = CONJ;
+  else {
+      PyErr_SetString(PyExc_ValueError, "trans must be N, T, or H");
+      return NULL;
+  }
+
+  if(L_N!=U_N){
+      PyErr_SetString(PyExc_TypeError, "L and U must have the same dimension");
+      return NULL;
+  }
+
+  if(!_CHECK_INTEGER(L_colind) || !_CHECK_INTEGER(L_rowptr) ||
+     !_CHECK_INTEGER(U_rowind) || !_CHECK_INTEGER(U_colptr) ){
+      PyErr_SetString(PyExc_TypeError, "row/column indices/pointers must be of type cint");
+      return NULL;
+  }
+
+  int L_type = PyArray_TYPE((PyArrayObject*)L_nzvals);
+  int U_type = PyArray_TYPE((PyArrayObject*)U_nzvals);
+  if (L_type != U_type) {
+      PyErr_SetString(PyExc_TypeError,
+                      "nzvals types of L and U differ");
+      return NULL;
+  }
+  if (!CHECK_SLU_TYPE(L_type)) {
+      PyErr_SetString(PyExc_TypeError,
+                      "nzvals is not of a type supported by SuperLU");
+      return NULL;
+  }
+
+  /* Create SuperLU matrices out of L and U. */
+  int* L_col_to_sup = intMalloc(L_N+1), *L_sup_to_col = intMalloc(L_N+1);
+  for(int i=0; i<=L_N; i++){
+    L_col_to_sup[i] = i;
+    L_sup_to_col[i] = i;
+  }
+  L_col_to_sup[L_N] = L_N-1;
+  SuperMatrix L_super = {0}, U_super={0};
+  int L_conv_err = SparseFormat_from_spMatrix(&L_super, L_N, L_N, L_nnz, -1,
+    (PyArrayObject*)L_nzvals, (PyArrayObject*)L_colind, (PyArrayObject*)L_rowptr, L_type, SLU_SC, SLU_TRLU, L_col_to_sup,L_sup_to_col);
+  if(L_conv_err){
+    return NULL;
+  }
+  int U_conv_err = SparseFormat_from_spMatrix(&U_super, U_N, U_N, U_nnz, 0,
+    (PyArrayObject*)U_nzvals, (PyArrayObject*)U_rowind, (PyArrayObject*)U_colptr, U_type, SLU_NC, SLU_TRU, NULL,NULL);
+  if(U_conv_err){
+    Destroy_SuperMatrix_Store((SuperMatrix*)&L_super);
+    return NULL;
+  }
+
+  /* Read right-hand-side/solution vector. */
+  PyArrayObject* X_arr = (PyArrayObject*)PyArray_FROMANY(
+      (PyObject*)X_py, L_type, 1, 2,
+      NPY_ARRAY_F_CONTIGUOUS | NPY_ARRAY_ENSURECOPY);
+  if (X_arr == NULL){
+      SUPERLU_FREE((void*)L_col_to_sup); SUPERLU_FREE((void*)L_sup_to_col);
+      Destroy_SuperMatrix_Store((SuperMatrix*)&L_super);
+      Destroy_SuperMatrix_Store((SuperMatrix*)&U_super);
+      return NULL;
+  }
+  if (PyArray_DIM((PyArrayObject*)X_arr, 0) != L_N) {
+      PyErr_SetString(PyExc_ValueError,
+                      "right hand side array has invalid shape");
+      SUPERLU_FREE((void*)L_col_to_sup); SUPERLU_FREE((void*)L_sup_to_col);
+      Destroy_SuperMatrix_Store((SuperMatrix*)&L_super);
+      Destroy_SuperMatrix_Store((SuperMatrix*)&U_super);
+      Py_DECREF(X_arr);
+      return NULL;
+  }
+
+  SuperMatrix X;
+  if (DenseSuper_from_Numeric((SuperMatrix*)&X, (PyObject*)X_arr)) {
+      SUPERLU_FREE((void*)L_col_to_sup); SUPERLU_FREE((void*)L_sup_to_col);
+      Destroy_SuperMatrix_Store((SuperMatrix*)&L_super);
+      Destroy_SuperMatrix_Store((SuperMatrix*)&U_super);
+      Py_DECREF(X_arr);
+      return NULL;
+  } /* X and X_arr share the same data but X_arr "owns" it. */
+
+  /* Call SuperLU functions. */
+  int info=0;
+  SuperLUStat_t stat = { 0 };
+  StatInit((SuperLUStat_t *)&stat);
+  int* perm_c = intMalloc(L_N);
+  for(int i=0; i<L_N; i++){
+    perm_c[i] = i;
+  }
+  int* perm_r = perm_c;
+
+  jmpbuf_ptr = (volatile jmp_buf *)superlu_python_jmpbuf();
+  SLU_BEGIN_THREADS;
+  if (setjmp(*(jmp_buf*)jmpbuf_ptr)) {
+      SLU_END_THREADS;
+      goto fail;
+  }
+  gstrs(L_type,
+        trans, &L_super, &U_super, perm_c, perm_r,
+        (SuperMatrix *)&X, (SuperLUStat_t *)&stat, (int *)&info);
+  SLU_END_THREADS;
+
+  if (info) {
+      PyErr_SetString(PyExc_SystemError,
+                      "gstrs was called with invalid arguments");
+      goto fail;
+  }
+
+  /* Deallocations and return. */
+  SUPERLU_FREE((void*)L_col_to_sup); SUPERLU_FREE((void*)L_sup_to_col);
+  SUPERLU_FREE((void*)perm_c);
+  Destroy_SuperMatrix_Store(&L_super);
+  Destroy_SuperMatrix_Store(&U_super);
+  XStatFree((SuperLUStat_t *)&stat);
+
+  return Py_BuildValue("Ni", X_arr, info);
+
+  fail:
+  SUPERLU_FREE((void*)L_col_to_sup); SUPERLU_FREE((void*)L_sup_to_col);
+  SUPERLU_FREE((void*)perm_c);
+  Destroy_SuperMatrix_Store(&L_super);
+  Destroy_SuperMatrix_Store(&U_super);
+  XStatFree((SuperLUStat_t *)&stat);
+  Py_DECREF(X_arr);
+  return NULL;
+}
+
 static char gssv_doc[] =
     "Direct inversion of sparse matrix.\n\nX = gssv(A,B) solves A*X = B for X.";
 
@@ -301,6 +466,15 @@ ilu : bool\n\
     (default: false)\n\
 ";
 
+static char gstrs_doc[] =
+  "gstrs(trans, L..., U..., b)\n"
+  "\n"
+  "solves the linear system A*x = b via backward and forward substitution, with\n"
+  "either A=L*U, A=(L*U)^T, or A=conj(L*U)^T.\n"
+  "trans says whether the matrix is transposed (\"T\"), conjugate transposed (\"H\"), or normal (\"N\").\n"
+  "L is specified in the compressed sparse column format via parameters N,nnz,nzvals,colind,rowptr.\n"
+  "U is specified in the compressed sparse column format via parameters N,nnz,nzvals,rowind,colptr.\n"
+  "b is specified as a dense vector.\n";
 
 /*
  * Main SuperLU module
@@ -311,6 +485,8 @@ static PyMethodDef SuperLU_Methods[] = {
      gssv_doc},
     {"gstrf", (PyCFunction) Py_gstrf, METH_VARARGS | METH_KEYWORDS,
      gstrf_doc},
+    {"gstrs", (PyCFunction) Py_gstrs, METH_VARARGS | METH_KEYWORDS,
+     gstrs_doc},
     {NULL, NULL}
 };
 
diff --git a/scipy/sparse/linalg/_dsolve/_superluobject.c b/scipy/sparse/linalg/_dsolve/_superluobject.c
index eaa125647560..d66375cfd74b 100644
--- a/scipy/sparse/linalg/_dsolve/_superluobject.c
+++ b/scipy/sparse/linalg/_dsolve/_superluobject.c
@@ -328,68 +328,37 @@ int NRFormat_from_spMatrix(SuperMatrix * A, int m, int n, int nnz,
                            PyArrayObject * nzvals, PyArrayObject * colind,
                            PyArrayObject * rowptr, int typenum)
 {
-    volatile int ok = 0;
-    volatile jmp_buf *jmpbuf_ptr;
-
-    ok = (PyArray_EquivTypenums(PyArray_TYPE(nzvals), typenum) &&
-          PyArray_EquivTypenums(PyArray_TYPE(colind), NPY_INT) &&
-          PyArray_EquivTypenums(PyArray_TYPE(rowptr), NPY_INT) &&
-          PyArray_NDIM(nzvals) == 1 &&
-          PyArray_NDIM(colind) == 1 &&
-          PyArray_NDIM(rowptr) == 1 &&
-          PyArray_IS_C_CONTIGUOUS(nzvals) &&
-          PyArray_IS_C_CONTIGUOUS(colind) &&
-          PyArray_IS_C_CONTIGUOUS(rowptr) &&
-          nnz <= PyArray_DIM(nzvals, 0) &&
-          nnz <= PyArray_DIM(colind, 0) &&
-          m+1 <= PyArray_DIM(rowptr, 0));
-    if (!ok) {
-        PyErr_SetString(PyExc_ValueError,
-                        "sparse matrix arrays must be 1-D C-contiguous and of proper "
-                        "sizes and types");
-        return -1;
-    }
-
-    jmpbuf_ptr = (volatile jmp_buf *)superlu_python_jmpbuf();
-    if (setjmp(*(jmp_buf*)jmpbuf_ptr)) {
-        return -1;
-    }
-    else {
-        if (!CHECK_SLU_TYPE(PyArray_TYPE(nzvals))) {
-            PyErr_SetString(PyExc_TypeError, "Invalid type for array.");
-            return -1;
-        }
-        Create_CompRow_Matrix(PyArray_TYPE(nzvals),
-                              A, m, n, nnz, PyArray_DATA((PyArrayObject*)nzvals),
-                              (int *) PyArray_DATA((PyArrayObject*)colind),
-                              (int *) PyArray_DATA((PyArrayObject*)rowptr),
-                              SLU_NR,
-                              NPY_TYPECODE_TO_SLU(PyArray_TYPE((PyArrayObject*)nzvals)),
-                              SLU_GE);
-    }
-
-    return 0;
+  return SparseFormat_from_spMatrix(A,m,n,nnz,1,nzvals,colind,rowptr,typenum,SLU_NR,SLU_GE,NULL,NULL);
 }
 
 int NCFormat_from_spMatrix(SuperMatrix * A, int m, int n, int nnz,
                            PyArrayObject * nzvals, PyArrayObject * rowind,
                            PyArrayObject * colptr, int typenum)
+{
+  return SparseFormat_from_spMatrix(A,m,n,nnz,0,nzvals,rowind,colptr,typenum,SLU_NC,SLU_GE,NULL,NULL);
+}
+
+int SparseFormat_from_spMatrix(SuperMatrix * A, int m, int n, int nnz, int csr /* 1=csr 0=csc -1=supernode csc */,
+                           PyArrayObject * nzvals, PyArrayObject * indices,
+                           PyArrayObject * pointers, int typenum, Stype_t stype, Mtype_t mtype,
+                           int* identity_col_to_sup, /* must be int array {0,1,2,...,n-2,n-1,n-1} */
+                           int* identity_sup_to_col /* must be int array {0,1,...,n} */)
 {
     volatile int ok = 0;
     volatile jmp_buf *jmpbuf_ptr;
 
     ok = (PyArray_EquivTypenums(PyArray_TYPE(nzvals), typenum) &&
-          PyArray_EquivTypenums(PyArray_TYPE(rowind), NPY_INT) &&
-          PyArray_EquivTypenums(PyArray_TYPE(colptr), NPY_INT) &&
+          PyArray_EquivTypenums(PyArray_TYPE(indices), NPY_INT) &&
+          PyArray_EquivTypenums(PyArray_TYPE(pointers), NPY_INT) &&
           PyArray_NDIM(nzvals) == 1 &&
-          PyArray_NDIM(rowind) == 1 &&
-          PyArray_NDIM(colptr) == 1 &&
+          PyArray_NDIM(indices) == 1 &&
+          PyArray_NDIM(pointers) == 1 &&
           PyArray_IS_C_CONTIGUOUS(nzvals) &&
-          PyArray_IS_C_CONTIGUOUS(rowind) &&
-          PyArray_IS_C_CONTIGUOUS(colptr) &&
+          PyArray_IS_C_CONTIGUOUS(indices) &&
+          PyArray_IS_C_CONTIGUOUS(pointers) &&
           nnz <= PyArray_DIM(nzvals, 0) &&
-          nnz <= PyArray_DIM(rowind, 0) &&
-          n+1 <= PyArray_DIM(colptr, 0));
+          nnz <= PyArray_DIM(indices, 0) &&
+          (csr ? (m+1) : (n+1)) <= PyArray_DIM(pointers, 0));
     if (!ok) {
         PyErr_SetString(PyExc_ValueError,
                         "sparse matrix arrays must be 1-D C-contiguous and of proper "
@@ -406,12 +375,29 @@ int NCFormat_from_spMatrix(SuperMatrix * A, int m, int n, int nnz,
             PyErr_SetString(PyExc_TypeError, "Invalid type for array.");
             return -1;
         }
-        Create_CompCol_Matrix(PyArray_TYPE(nzvals),
+        if(csr==1){
+          Create_CompRow_Matrix(PyArray_TYPE(nzvals),
                               A, m, n, nnz, PyArray_DATA(nzvals),
-                              (int *) PyArray_DATA(rowind), (int *) PyArray_DATA(colptr),
-                              SLU_NC,
+                              (int *) PyArray_DATA(indices), (int *) PyArray_DATA(pointers),
+                              stype,
                               NPY_TYPECODE_TO_SLU(PyArray_TYPE(nzvals)),
-                              SLU_GE);
+                              mtype);
+        } else if(csr==0) {
+          Create_CompCol_Matrix(PyArray_TYPE(nzvals),
+                                A, m, n, nnz, PyArray_DATA(nzvals),
+                                (int *) PyArray_DATA(indices), (int *) PyArray_DATA(pointers),
+                                stype,
+                                NPY_TYPECODE_TO_SLU(PyArray_TYPE(nzvals)),
+                                mtype);
+        } else if(csr==-1){
+          Create_SuperNode_Matrix(PyArray_TYPE(nzvals),
+                                 A, m, n, nnz, PyArray_DATA(nzvals),
+                                 (int *) PyArray_DATA(pointers), (int *) PyArray_DATA(indices),
+                                 (int *) PyArray_DATA(pointers), identity_col_to_sup, identity_sup_to_col,
+                                 stype,
+                                 NPY_TYPECODE_TO_SLU(PyArray_TYPE(nzvals)),
+                                 mtype);
+        }
     }
 
     return 0;
diff --git a/scipy/sparse/linalg/_dsolve/_superluobject.h b/scipy/sparse/linalg/_dsolve/_superluobject.h
index 2cf99db9568c..1f828d782dd4 100644
--- a/scipy/sparse/linalg/_dsolve/_superluobject.h
+++ b/scipy/sparse/linalg/_dsolve/_superluobject.h
@@ -56,6 +56,10 @@ int NRFormat_from_spMatrix(SuperMatrix *, int, int, int, PyArrayObject *,
 			   PyArrayObject *, PyArrayObject *, int);
 int NCFormat_from_spMatrix(SuperMatrix *, int, int, int, PyArrayObject *,
 			   PyArrayObject *, PyArrayObject *, int);
+int SparseFormat_from_spMatrix(SuperMatrix * A, int m, int n, int nnz, int csr /* 1=csr 0=csc */,
+                           PyArrayObject * nzvals, PyArrayObject * indices,
+                           PyArrayObject * pointers, int typenum, Stype_t stype, Mtype_t mtype,
+                           int* identity_col_to_sup, int* identity_sup_to_col);
 int LU_to_csc_matrix(SuperMatrix *L, SuperMatrix *U,
                      PyObject **L_csc, PyObject **U_csc,
                      PyObject *py_csc_construct_func);
@@ -153,6 +157,13 @@ jmp_buf *superlu_python_jmpbuf(void);
 #define Create_CompCol_Matrix_ARGS Create_CompRow_Matrix_ARGS
 #define Create_CompCol_Matrix_ARGS_REF Create_CompRow_Matrix_ARGS_REF
 
+
+#define Create_SuperNode_Matrix_ARGS \
+  SuperMatrix *a, int b, int c, int d, \
+  void* e, int* f, int* g, int* h, int* i, int* j, \
+  Stype_t k, Dtype_t l, Mtype_t m
+#define Create_SuperNode_Matrix_ARGS_REF a,b,c,d,e,f,g,h,i,j,k,l,m
+
 TYPE_GENERIC_FUNC(gstrf, void);
 TYPE_GENERIC_FUNC(gsitrf, void);
 TYPE_GENERIC_FUNC(gstrs, void);
@@ -160,5 +171,6 @@ TYPE_GENERIC_FUNC(gssv, void);
 TYPE_GENERIC_FUNC(Create_Dense_Matrix, void);
 TYPE_GENERIC_FUNC(Create_CompRow_Matrix, void);
 TYPE_GENERIC_FUNC(Create_CompCol_Matrix, void);
+TYPE_GENERIC_FUNC(Create_SuperNode_Matrix, void);
 
 #endif				/* __SUPERLU_OBJECT */
diff --git a/scipy/sparse/linalg/_dsolve/linsolve.py b/scipy/sparse/linalg/_dsolve/linsolve.py
index 454d153967ad..52edcf1cb04e 100644
--- a/scipy/sparse/linalg/_dsolve/linsolve.py
+++ b/scipy/sparse/linalg/_dsolve/linsolve.py
@@ -3,7 +3,7 @@
 import numpy as np
 from numpy import asarray
 from scipy.sparse import (issparse,
-                          SparseEfficiencyWarning, csc_matrix, csr_matrix)
+                          SparseEfficiencyWarning, csc_matrix, eye, diags)
 from scipy.sparse._sputils import is_pydata_spmatrix
 from scipy.linalg import LinAlgError
 import copy
@@ -617,8 +617,7 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
         Whether `A` is a lower or upper triangular matrix.
         Default is lower triangular matrix.
     overwrite_A : bool, optional
-        Allow changing `A`. The indices of `A` are going to be sorted and zero
-        entries are going to be removed.
+        Allow changing `A`.
         Enabling gives a performance gain. Default is False.
     overwrite_b : bool, optional
         Allow overwriting data in `b`.
@@ -626,8 +625,7 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
         If `overwrite_b` is True, it should be ensured that
         `b` has an appropriate dtype to be able to store the result.
     unit_diagonal : bool, optional
-        If True, diagonal elements of `a` are assumed to be 1 and will not be
-        referenced.
+        If True, diagonal elements of `a` are assumed to be 1.
 
         .. versionadded:: 1.4.0
 
@@ -651,9 +649,9 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
     Examples
     --------
     >>> import numpy as np
-    >>> from scipy.sparse import csr_matrix
+    >>> from scipy.sparse import csc_array
     >>> from scipy.sparse.linalg import spsolve_triangular
-    >>> A = csr_matrix([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
+    >>> A = csc_array([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
     >>> B = np.array([[2, 0], [-1, 0], [2, 0]], dtype=float)
     >>> x = spsolve_triangular(A, B)
     >>> np.allclose(A.dot(x), B)
@@ -661,20 +659,42 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
     """
 
     if is_pydata_spmatrix(A):
-        A = A.to_scipy_sparse().tocsr()
-
-    # Check the input for correct type and format.
-    if not (issparse(A) and A.format == "csr"):
-        warn('CSR matrix format is required. Converting to CSR matrix.',
+        A = A.to_scipy_sparse().tocsc()
+    
+    trans = "N"
+    if issparse(A) and A.format == "csr":
+        A = A.T
+        trans = "T"
+        lower = not lower
+
+    if not (issparse(A) and A.format == "csc"): 
+        warn('CSC or CSR matrix format is required. Converting to CSC matrix.',
              SparseEfficiencyWarning)
-        A = csr_matrix(A)
+        A = csc_matrix(A)
     elif not overwrite_A:
         A = A.copy()
 
-    if A.shape[0] != A.shape[1]:
+
+    M, N = A.shape
+    if M != N:
         raise ValueError(
             f'A must be a square matrix but its shape is {A.shape}.')
 
+    if unit_diagonal:
+        A = A.tolil()
+        A.setdiag(1)
+        A = A.tocsc()
+    else:
+        diag = A.diagonal()
+        if np.any(diag==0):
+            raise LinAlgError(
+                'A is singular: zero entry on diagonal.')
+        invdiag = 1/diag
+        if trans=="N":
+            A = A @ diags(invdiag)
+        else:
+            A = (A.T @ diags(invdiag)).T
+
     # sum duplicates for non-canonical format
     A.sum_duplicates()
 
@@ -683,61 +703,40 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
     if b.ndim not in [1, 2]:
         raise ValueError(
             f'b must have 1 or 2 dims but its shape is {b.shape}.')
-    if A.shape[0] != b.shape[0]:
+    if M != b.shape[0]:
         raise ValueError(
             'The size of the dimensions of A must be equal to '
             'the size of the first dimension of b but the shape of A is '
             '{} and the shape of b is {}.'.format(A.shape, b.shape))
 
-    # Init x as (a copy of) b.
-    x_dtype = np.result_type(A.data, b, np.float64)
-    if overwrite_b:
-        if np.can_cast(b.dtype, x_dtype, casting='same_kind'):
-            x = b
-        else:
-            raise ValueError(
-                'Cannot overwrite b (dtype {}) with result '
-                'of type {}.'.format(b.dtype, x_dtype))
-    else:
-        x = b.astype(x_dtype, copy=True)
+    A = A._asfptype()  # upcast to a floating point format
+    result_dtype = np.promote_types(A.dtype, b.dtype)
+    if A.dtype != result_dtype:
+        A = A.astype(result_dtype)
+    if b.dtype != result_dtype:
+        b = b.astype(result_dtype)
+    if not overwrite_b:
+        b = b.copy()
 
-    # Choose forward or backward order.
     if lower:
-        row_indices = range(len(b))
+        L = A
+        U = csc_matrix((N, N), dtype=result_dtype)
     else:
-        row_indices = range(len(b) - 1, -1, -1)
-
-    # Fill x iteratively.
-    for i in row_indices:
-
-        # Get indices for i-th row.
-        indptr_start = A.indptr[i]
-        indptr_stop = A.indptr[i + 1]
-
-        if lower:
-            A_diagonal_index_row_i = indptr_stop - 1
-            A_off_diagonal_indices_row_i = slice(indptr_start, indptr_stop - 1)
-        else:
-            A_diagonal_index_row_i = indptr_start
-            A_off_diagonal_indices_row_i = slice(indptr_start + 1, indptr_stop)
-
-        # Check regularity and triangularity of A.
-        if not unit_diagonal and (indptr_stop <= indptr_start
-                                  or A.indices[A_diagonal_index_row_i] < i):
-            raise LinAlgError(
-                f'A is singular: diagonal {i} is zero.')
-        if not unit_diagonal and A.indices[A_diagonal_index_row_i] > i:
-            raise LinAlgError(
-                'A is not triangular: A[{}, {}] is nonzero.'
-                ''.format(i, A.indices[A_diagonal_index_row_i]))
-
-        # Incorporate off-diagonal entries.
-        A_column_indices_in_row_i = A.indices[A_off_diagonal_indices_row_i]
-        A_values_in_row_i = A.data[A_off_diagonal_indices_row_i]
-        x[i] -= np.dot(x[A_column_indices_in_row_i].T, A_values_in_row_i)
-
-        # Compute i-th entry of x.
-        if not unit_diagonal:
-            x[i] /= A.data[A_diagonal_index_row_i]
+        L = eye(N, dtype=result_dtype, format='csc')
+        U = A
+        U.setdiag(0)
+
+    x, info = _superlu.gstrs(trans,
+                             N, L.nnz, L.data, L.indices, L.indptr,
+                             N, U.nnz, U.data, U.indices, U.indptr,
+                             b)
+    if info:
+        raise LinAlgError(
+            'A is singular.')
+
+    if not unit_diagonal:
+        invdiag = invdiag.reshape(-1, *([1]*(len(x.shape)-1)))
+        x = x * invdiag
 
     return x
+
diff --git a/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py b/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
index 6a7a4e45a32a..96fd034e3ce9 100644
--- a/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
+++ b/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
@@ -730,7 +730,7 @@ def test_zero_diagonal(self):
         rng = np.random.default_rng(43876432987)
         A = rng.standard_normal((n, n))
         b = np.arange(n)
-        A = scipy.sparse.tril(A, k=0, format='csr')
+        A = scipy.sparse.tril(A, k=0, format='csc')
 
         x = spsolve_triangular(A, b, unit_diagonal=True, lower=True)
 
@@ -741,12 +741,12 @@ def test_zero_diagonal(self):
         A = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0]], dtype=np.float64)
         b = np.array([1., 2., 3.])
         with suppress_warnings() as sup:
-            sup.filter(SparseEfficiencyWarning, "CSR matrix format is")
+            sup.filter(SparseEfficiencyWarning, "CSC or CSR matrix format is")
             spsolve_triangular(A, b, unit_diagonal=True)
 
     def test_singular(self):
         n = 5
-        A = csr_matrix((n, n))
+        A = csc_matrix((n, n))
         b = np.arange(n)
         for lower in (True, False):
             assert_raises(scipy.linalg.LinAlgError, spsolve_triangular, A, b, lower=lower)
@@ -774,29 +774,33 @@ def test_input_types(self):
     @pytest.mark.timeout(120)  # prerelease_deps_coverage_64bit_blas job
     @sup_sparse_efficiency
     def test_random(self):
-        def random_triangle_matrix(n, lower=True):
-            A = scipy.sparse.random(n, n, density=0.1, format='coo')
+        def random_triangle_matrix(n, lower=True, format="csr"):
+            A = scipy.sparse.random(n, n, density=0.1, format='lil')
             if lower:
-                A = scipy.sparse.tril(A)
+                A = scipy.sparse.tril(A,format="lil")
             else:
-                A = scipy.sparse.triu(A)
-            A = A.tocsr(copy=False)
+                A = scipy.sparse.triu(A,format="lil")
             for i in range(n):
                 A[i, i] = np.random.rand() + 1
+            if format=="csc":
+                A = A.tocsc(copy=False)
+            else:
+                A = A.tocsr(copy=False)
             return A
 
         np.random.seed(1234)
         for lower in (True, False):
-            for n in (10, 10**2, 10**3):
-                A = random_triangle_matrix(n, lower=lower)
-                for m in (1, 10):
-                    for b in (np.random.rand(n, m),
-                              np.random.randint(-9, 9, (n, m)),
-                              np.random.randint(-9, 9, (n, m)) +
-                              np.random.randint(-9, 9, (n, m)) * 1j):
-                        x = spsolve_triangular(A, b, lower=lower)
-                        assert_array_almost_equal(A.dot(x), b)
-                        x = spsolve_triangular(A, b, lower=lower,
-                                               unit_diagonal=True)
-                        A.setdiag(1)
-                        assert_array_almost_equal(A.dot(x), b)
+            for format in ["csr", "csc"]:
+                for n in (10, 10**2, 10**3):
+                    A = random_triangle_matrix(n, lower=lower)
+                    for m in (1, 10):
+                        for b in (np.random.rand(n, m),
+                                  np.random.randint(-9, 9, (n, m)),
+                                  np.random.randint(-9, 9, (n, m)) +
+                                  np.random.randint(-9, 9, (n, m)) * 1j):
+                            x = spsolve_triangular(A, b, lower=lower)
+                            assert_array_almost_equal(A.dot(x), b)
+                            x = spsolve_triangular(A, b, lower=lower,
+                                                   unit_diagonal=True)
+                            A.setdiag(1)
+                            assert_array_almost_equal(A.dot(x), b)
diff --git a/scipy/sparse/tests/test_array_api.py b/scipy/sparse/tests/test_array_api.py
index 88deec4554c1..50fb65f7f173 100644
--- a/scipy/sparse/tests/test_array_api.py
+++ b/scipy/sparse/tests/test_array_api.py
@@ -254,7 +254,7 @@ def test_spsolve(B):
 
 
 def test_spsolve_triangular():
-    X = scipy.sparse.csr_array([
+    X = scipy.sparse.csc_array([
         [1, 0, 0, 0],
         [2, 1, 0, 0],
         [3, 2, 1, 0],


Code review comments:
diff_hunk: @@ -774,29 +774,33 @@ def test_input_types(self):
     @pytest.mark.timeout(120)  # prerelease_deps_coverage_64bit_blas job
     @sup_sparse_efficiency
     def test_random(self):
-        def random_triangle_matrix(n, lower=True):
-            A = scipy.sparse.random(n, n, density=0.1, format='coo')
+        def random_triangle_matrix(n, lower=True, format="csr"):
+            A = scipy.sparse.random(n, n, density=0.1, format='lil')
path: scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
text: @user1:
This code can be accelerated using the new NumPy generators. As this gets executed several times, it is worth it in my opinion.

```python
import scipy
import numpy as np

rng = np.random.default_rng()

rvs = rng.random

n = 10**3

%timeit A = scipy.sparse.random(n, n, density=0.1, format='lil')  # 29.4 ms
%timeit A = scipy.sparse.random(n, n, density=0.1, format='lil', random_state=rng, data_rvs=rvs)  # 16.4 ms
```

@user1:
Also at least in `scipy.stats`, fixed random seeds are often preferred over new random data for every test run. Do we have general guidelines for this in scipy?

@author:
> This code can be accelerated using the new NumPy generators. 

Thanks, changed in 8aceb5ac9a26b77e966330560d1b2bee1f538774.

> Also at least in `scipy.stats`, fixed random seeds are often preferred over new random data for every test run. 

Also in [test_linsolve.py](https://github.com/scipy/scipy/blob/bf819ba58a6ccd51cc6d6dacf88ab18a80355148/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py), constant seeds are used in several places via `random.seed` or `random.RandomState`. I therefore chose to also fix the seed here, see 37ee66609dbbcc3852d7a1b1282e27cd187745d8. The particular random seed used here has been selected at random.

diff_hunk: @@ -774,29 +774,33 @@ def test_input_types(self):
     @pytest.mark.timeout(120)  # prerelease_deps_coverage_64bit_blas job
     @sup_sparse_efficiency
     def test_random(self):
-        def random_triangle_matrix(n, lower=True):
-            A = scipy.sparse.random(n, n, density=0.1, format='coo')
+        def random_triangle_matrix(n, lower=True, format="csr"):
+            A = scipy.sparse.random(n, n, density=0.1, format='lil')
             if lower:
-                A = scipy.sparse.tril(A)
+                A = scipy.sparse.tril(A,format="lil")
path: scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
text: @user1:
```suggestion
                A = scipy.sparse.tril(A, format="lil")
```

diff_hunk: @@ -651,30 +649,52 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
     Examples
     --------
     >>> import numpy as np
-    >>> from scipy.sparse import csr_matrix
+    >>> from scipy.sparse import csc_array
     >>> from scipy.sparse.linalg import spsolve_triangular
-    >>> A = csr_matrix([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
+    >>> A = csc_array([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
     >>> B = np.array([[2, 0], [-1, 0], [2, 0]], dtype=float)
     >>> x = spsolve_triangular(A, B)
     >>> np.allclose(A.dot(x), B)
     True
     """
 
     if is_pydata_spmatrix(A):
-        A = A.to_scipy_sparse().tocsr()
-
-    # Check the input for correct type and format.
-    if not (issparse(A) and A.format == "csr"):
-        warn('CSR matrix format is required. Converting to CSR matrix.',
+        A = A.to_scipy_sparse().tocsc()
+    
+    trans = "N"
+    if issparse(A) and A.format == "csr":
+        A = A.T
+        trans = "T"
+        lower = not lower
+
+    if not (issparse(A) and A.format == "csc"): 
+        warn('CSC or CSR matrix format is required. Converting to CSC matrix.',
              SparseEfficiencyWarning)
-        A = csr_matrix(A)
+        A = csc_matrix(A)
     elif not overwrite_A:
         A = A.copy()
 
-    if A.shape[0] != A.shape[1]:
+
+    M, N = A.shape
+    if M != N:
         raise ValueError(
             f'A must be a square matrix but its shape is {A.shape}.')
 
+    if unit_diagonal:
+        A = A.tolil()
+        A.setdiag(1)
+        A = A.tocsc()
+    else:
+        diag = A.diagonal()
+        if np.any(diag==0):
+            raise LinAlgError(
+                'A is singular: zero entry on diagonal.')
+        invdiag = 1/diag
+        if trans=="N":
path: scipy/sparse/linalg/_dsolve/linsolve.py
text: @user1:
```suggestion
        if trans == "N":
```

diff_hunk: @@ -683,61 +703,40 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
     if b.ndim not in [1, 2]:
         raise ValueError(
             f'b must have 1 or 2 dims but its shape is {b.shape}.')
-    if A.shape[0] != b.shape[0]:
+    if M != b.shape[0]:
         raise ValueError(
             'The size of the dimensions of A must be equal to '
             'the size of the first dimension of b but the shape of A is '
             '{} and the shape of b is {}.'.format(A.shape, b.shape))
 
-    # Init x as (a copy of) b.
-    x_dtype = np.result_type(A.data, b, np.float64)
-    if overwrite_b:
-        if np.can_cast(b.dtype, x_dtype, casting='same_kind'):
-            x = b
-        else:
-            raise ValueError(
-                'Cannot overwrite b (dtype {}) with result '
-                'of type {}.'.format(b.dtype, x_dtype))
-    else:
-        x = b.astype(x_dtype, copy=True)
+    A = A._asfptype()  # upcast to a floating point format
+    result_dtype = np.promote_types(A.dtype, b.dtype)
+    if A.dtype != result_dtype:
+        A = A.astype(result_dtype)
+    if b.dtype != result_dtype:
+        b = b.astype(result_dtype)
+    if not overwrite_b:
+        b = b.copy()
 
-    # Choose forward or backward order.
     if lower:
-        row_indices = range(len(b))
+        L = A
+        U = csc_matrix((N, N), dtype=result_dtype)
     else:
-        row_indices = range(len(b) - 1, -1, -1)
-
-    # Fill x iteratively.
-    for i in row_indices:
-
-        # Get indices for i-th row.
-        indptr_start = A.indptr[i]
-        indptr_stop = A.indptr[i + 1]
-
-        if lower:
-            A_diagonal_index_row_i = indptr_stop - 1
-            A_off_diagonal_indices_row_i = slice(indptr_start, indptr_stop - 1)
-        else:
-            A_diagonal_index_row_i = indptr_start
-            A_off_diagonal_indices_row_i = slice(indptr_start + 1, indptr_stop)
-
-        # Check regularity and triangularity of A.
-        if not unit_diagonal and (indptr_stop <= indptr_start
-                                  or A.indices[A_diagonal_index_row_i] < i):
-            raise LinAlgError(
-                f'A is singular: diagonal {i} is zero.')
-        if not unit_diagonal and A.indices[A_diagonal_index_row_i] > i:
-            raise LinAlgError(
-                'A is not triangular: A[{}, {}] is nonzero.'
-                ''.format(i, A.indices[A_diagonal_index_row_i]))
-
-        # Incorporate off-diagonal entries.
-        A_column_indices_in_row_i = A.indices[A_off_diagonal_indices_row_i]
-        A_values_in_row_i = A.data[A_off_diagonal_indices_row_i]
-        x[i] -= np.dot(x[A_column_indices_in_row_i].T, A_values_in_row_i)
-
-        # Compute i-th entry of x.
-        if not unit_diagonal:
-            x[i] /= A.data[A_diagonal_index_row_i]
+        L = eye(N, dtype=result_dtype, format='csc')
+        U = A
+        U.setdiag(0)
+
+    x, info = _superlu.gstrs(trans,
+                             N, L.nnz, L.data, L.indices, L.indptr,
+                             N, U.nnz, U.data, U.indices, U.indptr,
+                             b)
+    if info:
+        raise LinAlgError(
+            'A is singular.')
+
+    if not unit_diagonal:
+        invdiag = invdiag.reshape(-1, *([1]*(len(x.shape)-1)))
path: scipy/sparse/linalg/_dsolve/linsolve.py
text: @user1:
```suggestion
        invdiag = invdiag.reshape(-1, *([1] * (len(x.shape) - 1)))
```

diff_hunk: @@ -774,29 +774,33 @@ def test_input_types(self):
     @pytest.mark.timeout(120)  # prerelease_deps_coverage_64bit_blas job
     @sup_sparse_efficiency
     def test_random(self):
-        def random_triangle_matrix(n, lower=True):
-            A = scipy.sparse.random(n, n, density=0.1, format='coo')
+        def random_triangle_matrix(n, lower=True, format="csr"):
+            A = scipy.sparse.random(n, n, density=0.1, format='lil')
             if lower:
-                A = scipy.sparse.tril(A)
+                A = scipy.sparse.tril(A,format="lil")
             else:
-                A = scipy.sparse.triu(A)
-            A = A.tocsr(copy=False)
+                A = scipy.sparse.triu(A,format="lil")
             for i in range(n):
                 A[i, i] = np.random.rand() + 1
+            if format=="csc":
+                A = A.tocsc(copy=False)
+            else:
+                A = A.tocsr(copy=False)
             return A
 
         np.random.seed(1234)
         for lower in (True, False):
-            for n in (10, 10**2, 10**3):
-                A = random_triangle_matrix(n, lower=lower)
-                for m in (1, 10):
-                    for b in (np.random.rand(n, m),
-                              np.random.randint(-9, 9, (n, m)),
-                              np.random.randint(-9, 9, (n, m)) +
-                              np.random.randint(-9, 9, (n, m)) * 1j):
-                        x = spsolve_triangular(A, b, lower=lower)
-                        assert_array_almost_equal(A.dot(x), b)
-                        x = spsolve_triangular(A, b, lower=lower,
-                                               unit_diagonal=True)
-                        A.setdiag(1)
-                        assert_array_almost_equal(A.dot(x), b)
+            for format in ["csr", "csc"]:
+                for n in (10, 10**2, 10**3):
+                    A = random_triangle_matrix(n, lower=lower)
path: scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
text: @user1:
I am always a fan of pytest parametrization for such tests as it makes debugging much easier in my opinion. Could you adapt this code? I think it is very much historically grown.

@author:
Thanks for suggesting this! e645e4cc799141efe9a387a8f280f690f10ce462 makes the test parametric, removing the nested loops in a single testcase. Additionally, in 2f6802eacc6f357728be34ce7cd2e06f0be5d8ce, I've added tests for complex-valued matrices.

diff_hunk: @@ -651,30 +649,52 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
     Examples
     --------
     >>> import numpy as np
-    >>> from scipy.sparse import csr_matrix
+    >>> from scipy.sparse import csc_array
     >>> from scipy.sparse.linalg import spsolve_triangular
-    >>> A = csr_matrix([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
+    >>> A = csc_array([[3, 0, 0], [1, -1, 0], [2, 0, 1]], dtype=float)
     >>> B = np.array([[2, 0], [-1, 0], [2, 0]], dtype=float)
     >>> x = spsolve_triangular(A, B)
     >>> np.allclose(A.dot(x), B)
     True
     """
 
     if is_pydata_spmatrix(A):
-        A = A.to_scipy_sparse().tocsr()
-
-    # Check the input for correct type and format.
-    if not (issparse(A) and A.format == "csr"):
-        warn('CSR matrix format is required. Converting to CSR matrix.',
+        A = A.to_scipy_sparse().tocsc()
+    
+    trans = "N"
+    if issparse(A) and A.format == "csr":
+        A = A.T
+        trans = "T"
+        lower = not lower
+
+    if not (issparse(A) and A.format == "csc"): 
+        warn('CSC or CSR matrix format is required. Converting to CSC matrix.',
              SparseEfficiencyWarning)
-        A = csr_matrix(A)
+        A = csc_matrix(A)
     elif not overwrite_A:
         A = A.copy()
 
-    if A.shape[0] != A.shape[1]:
+
+    M, N = A.shape
+    if M != N:
         raise ValueError(
             f'A must be a square matrix but its shape is {A.shape}.')
 
+    if unit_diagonal:
+        A = A.tolil()
+        A.setdiag(1)
+        A = A.tocsc()
+    else:
+        diag = A.diagonal()
+        if np.any(diag==0):
path: scipy/sparse/linalg/_dsolve/linsolve.py
text: @user1:
```suggestion
        if np.any(diag == 0):
```

diff_hunk: @@ -774,29 +774,33 @@ def test_input_types(self):
     @pytest.mark.timeout(120)  # prerelease_deps_coverage_64bit_blas job
     @sup_sparse_efficiency
     def test_random(self):
-        def random_triangle_matrix(n, lower=True):
-            A = scipy.sparse.random(n, n, density=0.1, format='coo')
+        def random_triangle_matrix(n, lower=True, format="csr"):
+            A = scipy.sparse.random(n, n, density=0.1, format='lil')
             if lower:
-                A = scipy.sparse.tril(A)
+                A = scipy.sparse.tril(A,format="lil")
             else:
-                A = scipy.sparse.triu(A)
-            A = A.tocsr(copy=False)
+                A = scipy.sparse.triu(A,format="lil")
path: scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
text: @user1:
```suggestion
                A = scipy.sparse.triu(A, format="lil")
```