
[Sponsors] 
September 16, 2003, 21:38 
Block Tridiagonal and Block Bidiagonal solver

#1 
Guest
Posts: n/a

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. 

September 17, 2003, 20:01 
Re: Block Tridiagonal and Block Bidiagonal solver

#2 
Guest
Posts: n/a

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

August 2, 2011, 15:02 

#3 
Member
Join Date: Jul 2011
Location: US
Posts: 39
Rep Power: 7 
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.


August 2, 2011, 17:34 

#4  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 482
Rep Power: 13 
Quote:
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. 

August 3, 2011, 10:24 

#5 
Member
Join Date: Jul 2011
Location: US
Posts: 39
Rep Power: 7 
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.
__________________
CFD engineering resource 

August 3, 2011, 17:46 

#6  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 482
Rep Power: 13 
Quote:
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. Quote:
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. 

August 3, 2011, 18:36 

#7 
Member
Join Date: Jul 2011
Location: US
Posts: 39
Rep Power: 7 
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 "realworld" 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.
__________________
CFD engineering resource 

August 3, 2011, 19:04 

#8  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 482
Rep Power: 13 
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: Quote:


February 18, 2015, 17:53 

#9  
New Member
AJ
Join Date: Feb 2012
Location: Austin, TX
Posts: 10
Rep Power: 6 
Quote:


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Periodic Block tridiagonal solver  Abullah  Main CFD Forum  2  November 22, 2007 01:21 
block tridiagonal system solver for periodic bc  Abdullah  Main CFD Forum  8  February 23, 2007 14:15 
Block Tridiagonal Solver  Abdulhafid M. Elfaghi  Main CFD Forum  2  December 23, 2006 13:20 
Block tridiagonal matrix solver  Grace  Main CFD Forum  0  February 23, 2002 15:18 
Block Tridiagonal  John  Main CFD Forum  2  July 30, 2001 22:10 