[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
# License: SciPy-compatible

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.
>


More information about the Scipy-dev mailing list