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

[Benny Malengier]

2009/12/5 Simon Clift <ssclift at gmail.com>:
> On Wednesday 02 December 2009 15:48:53 David Warde-Farley wrote:
>> On 2-Dec-09, at 3:43 AM, Benny Malengier wrote:
>> > would be a type of sparse matrix one can manipulate. T
>> Note that converting to CSR and then adding elements would probably be
>> an inefficient way to do it.
>
> Yes, a tridiagonal matrix usually arises from 1D problems and finite
> differences for advection diffusion problems (which I need to solve).  Usually
> those kind of Newtonian problems are positive definite so you can LU factorize
> without pivoting.  My problem is particularly well conditioned, so I'm
> interested in that, done at top speed.

I do the same but with 3 components and a moving interface boundary.
So then it is banded + an entry in the last column for the interface.
So my remark was about  having code for the non-moving interface
boundary which would be banded (or with one component tridiagonal),
but when you move to a moving boundary (Landau transformation applied,
so the grid is moving with the internal grid fixed over the (0,1)
interval), you could take the code written for the banded system, and
just add the one column. To have such code, you would need to be able
to change from an implementation of a sparse banded matrix class to a
CSR sparse matrix, without the need to change the rest of the code.
Well, I did not work any details so perhaps it would not be that easy,
but that is the background of my request for having banded/tridiag as
type of sparse matrices.

> Banded usually arises from 2D,3D finite differences on regular grids.  That
> structure lends itself to multiplying out bands using vector routines.  That
> can be particularly efficient when you are building linear systems where the
> coefficients are changing or non-linear.

I don't follow here. If you have 2D nxn grid, then the bands are one
off the diagonal, and then n off the diagonal. The banded matrix
implementations don't allow for this as they are a fixed number of
lower and upper diagonals, so you need to use csr sparse matrix also
here, no? One could make a generic banded matrix implementation on top
of csr as I suggested, in which you say beforehand what specific off
diagonals can have values. That would be nice to do 2D/3D method of
lines.

By coincidence, I just have a mailing thread in the sundials user list
about such a thing:
http://old.nabble.com/Re%3A-CVODE-SUNDIALS-Feature-Wish-%3A-operator-splitting-td26645909.html
I obtained a private reply of a person who extended CVODE to work on
the sparse matrix class of the Meschach C library with the specific
goal of doing 2D/3D method of lines on top of CVODE.

Anyway, if you actually do convection-diffusion of one component, then
I would advise you to use CVODE, or via the pysundials package, or via
the higher level abstraction scikits.odes I wrote (am writing), I
don't think a tridiag solver will be able to beat an adaptive time
stepping BDF scheme. I have done that extensively, also with added
algebraic condition on mass conservation, with great results.

Benny

> CSR is most useful on irregular grids, as already noted.
>
> Each has got its use and changing formats, especially if the problems are
> large and time is of the essence, is usually a bad idea.
>
> --
> 1129 Ibbetson Lane
> Mississauga, Ontario
> L5C 1K9       Canada
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>