# [SciPy-dev] Sparse BLAS in scipy

Anand Patil anand.prabhakar.patil@gmail....
Tue Jan 15 12:51:11 CST 2008

```Hi all,

Sorry it took me so long to get to this, not even sure if it's needed
any longer, but I think the attached satisfies my earlier need for a
sparse triangular solver.

Here's my thinking: Gaussian elimination is trivial for sparse upper
triangular matrices, so I just wrote a wrapper for scipy.linsolve.splu
that makes sure splu gets an upper triangular matrix by transposing
the input matrix if necessary. Does this make sense?

Cheers,
Anand

--Begin--

# Author: Anand Patil, 2007

import scipy.sparse as sp
import scipy.linsolve as ls

def reverse_trans(trans):
if trans=='N':
return 'T'
else:
return 'N'

class sparse_trisolver(object):
"""__call__(b, trans='N' or 'T')"""

def __init__(self, chol_fac, uplo='U'):
self.uplo = uplo
if uplo=='U':
self.splu = ls.splu(chol_fac)
else:
self.splu = ls.splu(chol_fac.T)

def __call__(self, b, trans='N'):
if self.uplo=='L':
trans = reverse_trans(trans)
if len(b.shape)==1:
out = self.splu.solve(b, trans)
else:
out = b.copy()
for i in xrange(b.shape[1]):
out[:,i] = self.splu.solve(b[:,i], trans)
return out

# Test
if __name__=='__main__':
from numpy import asmatrix
from numpy.random import normal
from numpy.linalg import cholesky, solve

A_dense = asmatrix(normal(size=(5,5)))
A_dense = cholesky(A_dense.T*A_dense).T

A_csc = sp.csc_matrix(A_dense)
A_csr = sp.csr_matrix(A_dense)

U_csc = sparse_trisolver(A_csc)
U_csr = sparse_trisolver(A_csr)

L_csr = sparse_trisolver(A_csc.T, 'L')
L_csc = sparse_trisolver(A_csr.T, 'L')

B = normal(size=(5,30))

U_real = solve(A_dense, B)
L_real = solve(A_dense.T, B)

U_csc_sol = U_csc(B)
U_csr_sol = U_csr(B)
L_csc_sol = L_csc(B)
L_csr_sol = L_csr(B)

for upper in U_csc_sol, U_csr_sol:
print abs(upper - U_real).max()

for lower in L_csc_sol, L_csr_sol:
print abs(lower - L_real).max()

--End--

On Nov 28, 2007 7:01 PM, Anand Patil <anand.prabhakar.patil@gmail.com> wrote:
> On Nov 28, 2007 5:55 PM, Anand Patil <anand.prabhakar.patil@gmail.com>
> wrote:
>
> >
> >
> >
> > OK, if you've done the research already and come to that conclusion I
> should just do it! Regardless of the backend (LU, gauss-seidel or from
> scratch) here's what I'm proposing for the Python interface:
> >
> > B=trimult(A, B, uplo='U', side='L', inplace=True)
> > B=trisolve(A, B, uplo, side, inplace)
> >
>
> I just sat down to write the 'trimult' routines and realized there's no
> need, the sparse general matrix-matrix multiplication saves the work anyway!
> Sorry for being slow. I'll just work on the 'trisolve' one.
>
```