[SciPy-user] Tridiagonal solver extensive use

[A. M. Archibald]

On 09/11/06, Ganesh V <ganesh.iitm at gmail.com> wrote:
> Hi!
>
>   I have a specific problem (actually fluid mechanics.. solving 2 species
> and 1 energy eqn for those concerned !!!) where I'll have to solve around 8
> eqns of the type
>
>                Ax =b where A is tridiagonal
>
>  PER TIME STEP to get the various derivatives.. And I may have to solve for
> a million time steps or so (pretty small.. of the order of microseconds). I
> was planning to use Scipy.lianalg.solve_banded for this purpose. How do I
> optimize storing these arrays if all/most of the diagonal/off diagonal terms
> are the same. I guess this comes in many places esp.. wherever Finite
> differencing is used. In my case, most of the Diagonal terms are one (except
> the first 2 and last 2) and same is the case with the off diagonal terms as
> well.
>
> Now comes the main question. Do I really go for doing this in Fortran or
> make my life easier by writing it in Python- Scipy. Actually I may finish
> writing the code pretty soon if it were in Python.. some 50-60 lines..!! ,
> but my only worry is performance !! Once the code is written I'll want
> results fast (I guess everybody does !!). Will I really gain by coding in
> Fortran (actually my programming skills aren't that great, so that may
> offset the benefit as well) ??

I would suggest implementing a tridiagonal solver directly in python
(it should be simple) if you can do so without using a for loop. It
probably won't really take advantage of the simplicity of your
coefficients, but it will run fairly quickly. If you have a solution
to your problem that works but is a bit too slow, then you can try
weave, pyrex, or writing a module in C or FORTRAN. That will certainly
be faster, but you don't necessarily care.

The general rule is, get it working first, then make it fast.

A. M. Archibald