[SciPy-dev] Tri-diagonal LAPACK Routines - Shall I interface them?

Benny Malengier benny.malengier@gmail....
Tue Dec 8 07:55:25 CST 2009


2009/12/7 Nathan Bell <wnbell@gmail.com>:
> On Wed, Dec 2, 2009 at 3:43 AM, Benny Malengier
> <benny.malengier@gmail.com> wrote:
>>
>> Interesting, this is exactly the function I needed for my problem, but
>> I was looking in scipy.sparse.linalg, so did not notice banded matrix
>> solver was present in scipy.linalg.
>>
>> In my logic, the "matrix diagonal orded form" of
>> http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html#scipy.linalg.solve_banded
>> would be a type of sparse matrix one can manipulate. This would allow
>> things like changing matrix diagonal orded form sparse matrix to a csr
>> matrix, adding some extra elements off the diagonals, and then calling
>> a more generic solver.
>>
>
> Can't you do that already with scipy.sparse.dia_matrix?  If I'm not
> mistaken, dia_matrix is (slightly) more general than the banded format
> but similarly efficient.
>
> In an ideal world scipy.sparse.spsolve() would detect the case that A
> was a dia_matrix (with small bandwidth) and invoke the LAPACK method
> in scipy.linalg instead of using the general sparse LU solver.

Yes, dia can do that, but Lapack expects another dia format, so there
would be some undesired overhead (see
http://www.netlib.org/lapack/double/dgbsv.f or
http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html#scipy.linalg.solve_banded)

On the other hand, umfpack desires csr, so from that point of view,
off diagonal dia matrices as stored by dia, would have to be converted
in the csr format as I suggested to implement a dia sparse matrix
directly. The present dia implementation would do A.to_coo().to_csr()
before a solution can be obtained from umfpack.

Benny

>
> --
> Nathan Bell wnbell@gmail.com
> http://www.wnbell.com/
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>


More information about the SciPy-Dev mailing list