Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 6 down vote favorite
1

I'm trying to solve a 2D poisson equation by finite differences. In the process, I obtain a sparse matrix with only 5 variables in each equation, For example, if the variable were U, then the discretization would yield:

$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$

I know that I can solve this system by an iterative method, but the thought occurred to me that if I ordered the variables appropriately, I might be able to obtain a banded matrix which could be solved by a direct method (i.e. gaussian elimination w/o pivoting). Is this possible? Are there any strategies for doing this for other, perhaps less structured sparse systems?

share|improve this question
1  
Something like Cuthill-McKee, then? – J. M. Jan 18 '12 at 17:46
Interesting... i've never heard of the Cuthill-McKee algorithm before! :) – Paul♦ Jan 18 '12 at 18:18
1  
There's also a Reverse Cuthill-McKee as well. – Geoff Oxberry♦ Jan 19 '12 at 4:38
1  
I hope it is clear from the answers, but you do not want to use a banded solver for this problem, nor choose an ordering that minimizes bandwidth. Perhaps the question or the chosen answer can be edited to make this clear, otherwise I fear that this myth will be perpetuated. I gave a visual comparison and compared fill in scicomp.stackexchange.com/a/880/119. – Jed Brown Jan 19 '12 at 13:40
@JedBrown: Actually, I'm not quite working with a poisson problem, per se... My problem has a similar structure to the poisson problem... The indicies of the variables (i's and j's) are exactly the same, and the matrix is diagonally dominant with the off-diagonal entries (within the same row) add to exactly the sum of the diagonal entry. – Paul♦ Jan 19 '12 at 14:09
show 1 more comment

4 Answers

up vote 10 down vote accepted

This is a well-studied problem in the field of sparse-direct solvers. I highly recommend reading Joseph Liu's overview of the multifrontal method in order to get a better idea of how reorderings and supernodes effect fill-in and solution time.

Nested dissection is an extremely common way to generate the reordering, and essentially consists of recursive graph partitioning. MeTiS is the de facto standard for graph partitioning, and you can read about some of the ideas behind it here. Another commonly used package is SCOTCH, and Chaco is also important, as its authors introduced multi-level graph partitioning, which is also the fundamental idea behind MeTiS.

George and Liu showed in their classic book that 2d sparse-direct solutions only require $O(n^{3/2})$ work and $O(n \log n)$ memory, while 3d sparse-direct requires $O(n^2)$ work and $O(n^{4/3})$ memory.

share|improve this answer
Do you have a citation for the George and Liu reference? – Paul♦ Jan 18 '12 at 18:49
Added; I was about to get out of the car when I first submitted it. I know that there exists a freely available version of the book online somewhere (Jed knows where it is), but I could not find it. – Jack Poulson Jan 18 '12 at 19:04
I updated the link to point to the PDF of the book instead of the book review. – Jed Brown Jan 18 '12 at 20:56
@JedBrown That was a great reference! Thanks so much! :) – Paul♦ Jan 19 '12 at 5:34
1  
@Alexander Everyone attributes the 3D bound to George and Liu, though I don't know if they explicitly point it out in the book. It is obvious from the theory however. The minimal vertex separator for a $n = m\times m\times m$ grid is $n^{2/3} = m\times m$. The dense matrix associated with that supernode has $(n^{2/3})^2 = n^{4/3}$ entries and requires $(n^{2/3})^3 = n^2$ operations to factor. The logarithmic term in the 2D case is more subtle and is treated in Chapter 8 on Nested Dissection, which achieves the lower bound. – Jed Brown Aug 27 '12 at 13:59
show 1 more comment
up vote 4 down vote

Speaking of multifrontal methods, Tim Davis, who works on multifrontal methods for LU factorization (UMFPACK) has a number of routines that will reorder matrices to minimize fill-in. You can find them as here as part of SuiteSparse. SuiteSparse uses MeTiS.

One other thing to note: In some problems, you can be clever about ordering variables so that you get banded, or close to banded, patterns, which can save you the trouble (and the CPU time) of calling these algorithms. However, this clever reordering requires insight on your part and is nowhere near as general as the graph-theory-based reordering algorithms people have mentioned in their answers here.

share|improve this answer
Awesome Geoff! Thanks! :) – Paul♦ Jan 19 '12 at 4:52
You're welcome, Paul. If you like it, vote it up. – Geoff Oxberry♦ Jan 19 '12 at 5:31
up vote 3 down vote

There's an algorithm called ADI (Alternating Direction Implicit) in applied math circles and Split-operator in physics circles that does basically what you describe. It's an iterative method, and it follows this basic procedure:

  1. For every value of $y$ , relax in the $x$-direction. This matrix should be tridiagonal, so it can be solved directly in relatively little time.

  2. For every value of $x$ , relax in the $y$-direction. Again, this should be pretty quick.

  3. Repeat 1 and 2 until the error is as small as you want it to be.

I don't know the formal complexity of this algorithm, but I've found it to converge in fewer iterations than things like Jacobi and Gauss-Seidel every time I've used it.

share|improve this answer
2  
If you decide to go the operator splitting route, something you'll want to be careful about is that operator splitting is known to result in errors in steady state solutions in some cases. (One of my labmates has developed a way to overcome this difficulty, but I don't believe he's published it yet.) Also, operator splitting is known to cause numerical errors. There are well-established ways to estimate these errors a posteriori; Don Estep has done excellent work in that area. – Geoff Oxberry♦ Jan 19 '12 at 4:37
@GeoffOxberry It sounds like you are referring to a different splitting. You can use ADI in a fully implicit scheme which has no splitting error because it actually solves the system. There are also IMEX methods that rigorously control splitting errors. – Jed Brown Jan 19 '12 at 5:14
@JedBrown: I was talking about Godunov and Strang splitting, which can yield similar tridiagonal matrices if you split up the $x$ and $y$ terms. My bad. (I keep telling people I want to learn PDE methods in a postdoc...) – Geoff Oxberry♦ Jan 19 '12 at 5:24
I've never heard of Godunov and Strang splitting. I tend to split my operator with the Baker-Campbell-Hausdorf formula. Is that the same thing? – Dan Jan 19 '12 at 6:20
up vote 2 down vote

Cuthill-Mckee is the de-factor standard for what you want to do. If you wanted to play with this method, there's an easy to use implementation of the algorithm (and its reverse) in the BOOST Graph library, and the documentation contains examples how to use it.

share|improve this answer
actually reverse Cuhill-McKee; it usually gives less fill. But a nested dissection ordering is far superior to a low bandwidth ordering. – Arnold Neumaier Mar 29 '12 at 19:46

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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