I am solving a physical problem using implicit numerical scheme. This leads me to solving a linear equation with tridiagonal matrix. I've coded this algorithm from Wikipedia. I wonder if there is an efficient library which allows to solve this type of equation in optimised way. An important note is that matrix itself changes only when system parameters are changing, so I had an opportunity to precalculate some algorithm steps for a nice peformance bonus. I am using C++.
$\begingroup$
$\endgroup$
2
-
$\begingroup$ How big of a system is it, does it need to be parallel? $\endgroup$– aterrelCommented Nov 30, 2011 at 22:59
-
1$\begingroup$ Size depends on required accuracy (from hundred to tens of thousands values). Now I am coding on one-core computer, but it's possible to get access to university supercomputer with many cpus available, so parallelism support would be nice. $\endgroup$– gmkCommented Nov 30, 2011 at 23:17
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
2
You should probably start with the LAPACK implementation, ?gtsv, e.g., dgtsv. If you want a distributed-memory version, then you might want to start with ScaLAPACK's p?gtsv.
EDIT: Since your matrix does not change very often, you can avoid redundantly factoring the tridiagonal matrix by breaking up the LAPACK routine ?gtsv into the factorization step, ?gttrf, and the solve stage, ?gttrs. Similarly named routines exist in ScaLAPACK that serve the same purpose.
-
$\begingroup$ Thank you, it looks like what I need. I'll try now to run this routines from my code. $\endgroup$– gmkCommented Nov 30, 2011 at 23:59
-
1$\begingroup$ Since you are calling it from C++, make sure to declare the prototype inside of an extern "C" { } block. Depending on your system, you may need to append an underscore to the routine name. $\endgroup$ Commented Dec 1, 2011 at 0:06
$\begingroup$
$\endgroup$
For distributed parallel systems: I have not tried ScaLAPACK, which has a parallel tridiagonal solver, for which there are examples available online. I have tried with some success a method proposed by David Moulton in a LANL publication. Coding this up might be more than you want to do, but by using LAPACK, it is strait forward.
$\begingroup$
$\endgroup$
1
I found an interesting recursive algorithm here at page 975. It looks promising, I wonder what more experienced people say about it.
-
$\begingroup$ Numerical Recipes has some errors in it. In terms of a source of codes to use, it isn't the best, although some consider it a classic. I'd be surprised if ScaLAPACK didn't implement an algorithm at least as efficient as recursive cyclic reduction. $\endgroup$ Commented Mar 6, 2012 at 10:49