Hi,

Can anybody tell me where I can find an effecient Block Tridiagonal and a Block Bidiagonal Solver in C. I need this to solve implicit Euler equations. Thanks in advance. Amit.

Have you tried the Numerical Recipes?... Look for it at the IMSL Fortran Library (Fortran 90)

Hope that helps you.

[Docfreezzzz]

There is a Block Tridiagonal Solver available with full source code here that I have used for many CFD projects. This is mostly applicable to 1D CFD problems, however.

[arjun]

Originally Posted by Docfreezzzz

There is a Block Tridiagonal Solver available with full source code here that I have used for many CFD projects. This is mostly applicable to 1D CFD problems, however.

This is based on Thomas' algorithm so in case of block matrices it is NOT efficient. If your case is large lets say block of 2000 * 2000 , it will take lot of time.

If possible one should try Block Cyclic reduction methods, they are very efficient for block matrices.

[Docfreezzzz]

You are correct. However, that method is much more complex and generally harder to follow for the inexperienced. Also, may I suggest that for most tridiagonal problems the additional complexity is not warranted for such a small time savings. Any code that I have ever needed the capability to solve a block tridiagonal system is so fast to begin with that an increase in speed would be almost imperceptible.

[arjun]

Originally Posted by Docfreezzzz

You are correct. However, that method is much more complex and generally harder to follow for the inexperienced. Also, may I suggest that for most tridiagonal problems the additional complexity is not warranted for such a small time savings. Any code that I have ever needed the capability to solve a block tridiagonal system is so fast to begin with that an increase in speed would be almost imperceptible.

If take say mesh of 1000 x 1000 ( or say n x m )and discretesize Poisson problem on it.

The operation count in block cyclic reduction is roughly nmlog(m) .

But if you use the method of Thomas, you are looking at inverting m x m system n times . And have to do it 2 times because of forward substitution and then backward substitution. Now with LU factorization you operation count for inverting the system would be say 3m^2 So total cost is

(3m^2 ) x n x 2 + other operations like matrix matrix multiplications etc.

For 1000x 1000 mesh you are looking at

3 x 1000 x 1000 x 1000 x 2

compared to 1000 x 1000 x log2(1000) of block cyclic reduction.

The ratio by which block cyclic reduction is faster (roughly) : 6000 / log2(1000) .

Which I would not call close. In fact many iterative methods would be much faster than what Thomas algorithm would achieve in this case.

Originally Posted by Docfreezzzz

You are correct. However, that method is much more complex and generally harder to follow for the inexperienced.

Yes, writing that code is difficult, so easy way out is to download fishpack and use the routine from there.
It is really difficult to code it yourself. I did it but would not recommend anyone to do it.





PS: There are approaximate block cyclic reduction algorithms exists, they ,instead of
nmlog2(m) , can do in nklog2(m) where k ~ 10 or 15.

That means for uniform mesh I could reduce the timing from 1000x1000xlog2(1000) to 15x1000xlog2(1000) . All this is different story though.

[Docfreezzzz]

The cost of solving the discrete Euler equations on a 1000 node 1D grid (which results in a block tridiagonal system) with the algorithm I posted is more like 3*1000*(3^2), since the block sizes are only 3x3. This makes the situation far less dire than you propose with the Poisson equation. However, if you read my post carefully I was not disagreeing with you in the slightest. I was merely stating what a feel is a near universal truth. "Premature optimization is pure evil". That is, if you don't know what the costs are to begin with, in terms of wall clock time, then optimizing a routine based on theoretical operations counts will only get you into trouble.

Often this sort of activity results in an algorithm that is only moderately faster in practice due to cache trashing, etc. and a whole boatload of undecipherable code for the next in line to deal with. Also, most 1D Euler codes which is what I assume the OP is writing are often not of much interest past a few thousand nodes/cvs. In this case, I would propose that the additional effort is simply not worth the end result.

As for iterative methods, these are the only shot most of us get when dealing with "real-world" CFD problems due to the expense of direct methods. In fact, were the OP developing a more complex code I would suggest dropping direct methods altogether as they are completely unnecessary and often unusable on large problems due to memory limitations.

Kudos. Thanks for the education on block cyclic methods. I had never had a need for these but I'll keep that in my toolbox should it ever arise.

[arjun]

I give you that inverting 3x3 matrix is not very expensive. But it seems you missed my point. Which was Thomas algorithm is not very efficient for general block tridiagonal system. It could be used for smaller blocks where inverting them is not an issue.

This is why I said in my first post:

Originally Posted by arjun

This is based on Thomas' algorithm so in case of block matrices it is NOT efficient. If your case is large lets say block of 2000 * 2000 , it will take lot of time.

If possible one should try Block Cyclic reduction methods, they are very efficient for block matrices.

[akhi2908]

Originally Posted by Docfreezzzz

There is a Block Tridiagonal Solver available with full source code here that I have used for many CFD projects. This is mostly applicable to 1D CFD problems, however.

I want to use this code. Can you please update the link?