solve_triangular vs spsolve_triangular: Treating non-triangular matrices

[nschloe]

I need to solve upper-triangular systems with both dense and sparse matrices. To this end, I'm looking at solve_triangular for dense and spsolve_triangular for sparse.

solve_triangular correctly discards the unused part, i.e., the solution works even if the input matrix A is not triangular. This is useful in many cases, e.g., in stationary methods like Gauss-Seidel; it saves you all matrix-manipulation shenanigans.

import numpy as np
from scipy.linalg import solve_triangular

A = np.array([
    [1.0, 2.0],
    [3.0, 4.0]
])
b = np.array([1.0, 1.0])

sol = solve_triangular(A, b)
print(sol)

[0.5, 0.25]

spsolve_triangular one the other hand flat out won't accept non-triangular matrices:

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve_triangular


A = csr_matrix([
    [1.0, 2.0],
    [3.0, 4.0]
])
b = np.array([1.0, 1.0])

sol = spsolve_triangular(A, b)

numpy.linalg.LinAlgError: A is not triangular: A[0, 1] is nonzero.

This is a feature request for accepting non-triangular matrices to spsolve_triangular to avoid having to create temporary matrices for the solution of triangular problems.

This is with SciPy 1.6.3.

[perimosocordiae]

Here's where we do the check:

scipy/scipy/sparse/linalg/dsolve/linsolve.py

Line 613 in 701ffcc

if A.indices[A_diagonal_index_row_i] > i:

It's a little complex because we're doing indexing tricks to avoid seeking through the nonzero entries of each row to find the diagonal entry. I think a rewrite could improve the code's readability and remove this limitation at the same time. We'll want to be careful not to incur a performance penalty, as well.