9
$\begingroup$

In my project I have to solve a couple of tridiagonal matrices at every time step, so it is crucial to have a good solver for those. I did my own implementation, just the classical way to do it described on Wikipedia. I then tried using Lapack instead, and to my surprise it was slower!

Now, inside Lapack it seems like it does the solving by LU factorization and I wonder why, isn't it more complex than it could be?

Additionally, I found an algorithm in the "Numerical Recipes" book from nr.com which recursively divides the system into smaller tridiagonal problems. It looked promising. Are there any other goodies out there?

Update: the problem size is about 1000x1000. I used GotoBLAS, it gives you a Lapack 3.1.1 library also. The problem is not symmetric. I used the Lapack routine for general tridiagonal matrices.

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ You will have to state which LAPACK routines you used for this. Note that dgtsv performs partial pivoting, but your code may not do this. Please also state which LAPACK implementation you tested with and what problem sizes you benchmarked with. Also, is your problem symmetric positive definite? $\endgroup$
    – Jed Brown
    Commented Mar 5, 2012 at 15:53
  • $\begingroup$ I added some info in the question formulation. $\endgroup$
    – tiam
    Commented Mar 5, 2012 at 16:01
  • $\begingroup$ Is your application something to do with Finite Volume Methods? $\endgroup$
    – Inquest
    Commented Mar 5, 2012 at 16:32
  • $\begingroup$ It is finite differences, but in this perspective it is more or less the same I guess. $\endgroup$
    – tiam
    Commented Mar 5, 2012 at 16:54

1 Answer 1

Reset to default
7
$\begingroup$

You are using a reference implementation that does partial pivoting. Tridiagonal solves do very little work and do not call into the BLAS. It is likely slower than your code because it does partial pivoting. The source code for dgtsv is straightforward.

If you will solve with the same matrix multiple times, you may want to store the factors by using dgttrf and dgttrs. It is possible that the implementations in an optimized LAPACK such as MKL, ACML, or ESSL will be more performant.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I'm a bit curious. Gaussian Elim with PP would work for all matrices including TriDiagonal. In CFD, we use a special method for FVM 1D cases called TDMA. Which do you reckon would be faster for the case he is discussing? Although, I am not entirely sure his matrices are diagonally dominant. $\endgroup$
    – Inquest
    Commented Mar 5, 2012 at 16:31
  • $\begingroup$ The TDMA is what I implemented in my code. The question is why would the super-fast Lapack use the partial pivoting procedure in such a particular matrix, which is solved faster by such an easy method as TDMA. $\endgroup$
    – tiam
    Commented Mar 5, 2012 at 16:44
  • 1
    $\begingroup$ It is exactly the same algorithm (Gaussian elimination specialized for a tridiagonal matrix), but your implementation does not do partial pivoting, so it may be numerically unstable. That pivoting is not free and your are comparing to the reference implementation. The reference implementation is not optimized for performance and the partial pivoting is not free. $\endgroup$
    – Jed Brown
    Commented Mar 5, 2012 at 17:11
  • $\begingroup$ I see what you mean, I get advantage from my knowledge about the systems I am solving. Do other implementations of LAPACK give performance boosts due to adaptation to specific architecture or does it go beyond that? $\endgroup$
    – tiam
    Commented Mar 5, 2012 at 17:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.