For some time now I've been missing a function in scipy that exploits the triangular structure of a matrix to efficiently solve the associated system, so I decided to implement it by binding the LAPACK method "trtrs", which also checks for singularities and is capable handling several right-hand sides. Contrary to what I expected, binding Fortran code with f2py is pretty straightforward, even for someone like me who has never programmed in that language: I took a similar example, modified it's parameters and it worked! Also, thanks to Pauli Virtanen the review process was really fast and the patch was committed within a few hours. The high level interface for LAPACK's trtrs is linalg.solve_triangular, which accepts roughly the same arguments as linalg.solve, but assumes the first argument is a triangular matrix: [cc lang="python"] In [1]: from scipy import linalg In [2]: linalg.solve_triangular([[1, 1], [0, 1]], [0, 1]) Out[2]: array([-1., 1.]) [/cc] Simple benchmarks lets us clearly appreciate the complexity gap between both methods : solving an (n, n) triangular system is an O(n^2) operation, while solving a full one is at least a O(n^3):