2

Since the computational complexity of direct elimilation methods for solving linear systems is $O(n^3)$, it's not practical when the number of dofs is large. But how large would you call it a large system? I saw some material saying that 100,000 dofs would be the watershed when you should consider iterative methods. Wouldn't that be too large? I mean when dofs reach just less than 100,000, the complexity would be $O(100,000^3)$ with direct methods, it's hard to imagine that we can afford such a huge comlexity? Could anyone give some concrete examples or intuitive sugestions indicating the usual upper bound when direct elimilation methods are treated as impractical? Literatures would also be fine. Thank you!

user123
  • 679
  • 1
  • 5
  • 16
  • 6
    It's impossible to give such a bound, since the actual (not theoretical worst-case!) performance depends a) on the structure of the matrix (extreme example: triangular) and b) the hardware (which is a moving target). Also, for direct solvers, the available memory is in general the bottleneck, so the general rule of thumb is that if the factors fit into your memory, you can use a direct solver, and if not (and you usually don't know until you try), use an iterative solver. – Christian Clason Nov 20 '16 at 10:26
  • Thank you @Christian, but what if I restrict the matrix to those resulting from FEM, so it would be sparse, but not necessarily banded due to ordering issues. Then what would be your considerations? – user123 Nov 20 '16 at 10:51
  • The same, since the structure of the matrix still very much depends on the PDE you're trying to solve and the point about the hardware still stands -- the answer will be very different for a PC, an HPC cluster, or a GPU. – Christian Clason Nov 20 '16 at 10:54
  • OK, let's say it's on a personal computer and we are trying to solve the elliptic PDE. At least the $O(n^3)$ could give some guidance? – user123 Nov 20 '16 at 11:06
  • 4
    No, since that's only a theoretical worst-case estimate and doesn't tell you anything about a specific matrix. My guidance in your situation would be to try a direct solver, and if you get an "out of memory" error, switch to an iterative solver. – Christian Clason Nov 20 '16 at 11:08
  • In finite element structural analysis, direct solvers are almost always used as opposed to iterative solvers. Solving a "typical" 100000-equation FE model on a basic PC is routine. Slide two in this presentation from ABAQUS defines a "large" FE model as one with $> 5M$ equations (http://imechanica.org/files/l10-large-models_2.pdf). Note that this presentation is from 2005 so that estimate is surely larger now. – Bill Greene Nov 20 '16 at 14:02
  • So you mean for these FE model with > 5M equations, it is prefered to use the direct solvers? From my knowledge, direct solvers for these large sparse system would introduce some fill-in issues and thus not quite effective as iterative solver, although some reoredering techniques may be utilized. – user123 Nov 20 '16 at 14:11
  • 2
    @David That's the point: You can't generalize -- every problem is different, and there's a huge variety of both direct and iterative solvers (not to mention preconditioners) that behave differently for different problems (and different hardware). For some of the issues, see this question: http://scicomp.stackexchange.com/questions/3262/how-to-choose-a-method-for-solving-linear-equations – Christian Clason Nov 20 '16 at 14:43
  • 2
    As a concrete example, I just now solved the equations from a $70 \times 70 \times 70$ discretization (343K equations) of the Laplace equation using a direct solver (MUMPS) on a 6-year-old, moderate-power PC with 12 GB of memory in approximately one minute. As an aside, I have been puzzled for many years about the misconceptions regarding state-of-the-art sparse direct solvers. – Bill Greene Nov 20 '16 at 14:54
  • Wow, that's impressive. To be honest, I don't have much experience with these solvers in such large practical problems and that's also why I have this question. – user123 Nov 20 '16 at 14:58
  • 1
    With good orderings, the complexity of factorization generally scales with the cost of doing dense operations on the "separator fronts" -- collections of degrees of freedom that, when removed, partition the remaining degrees of freedom into two noninteracting pieces (and so on, recursively). So, doing a sparse factorization of a 2D problem is, up a constant and some logarithmic factor, like doing a dense factorization of a 1D problem, and doing a sparse factorization of a 3D problem is like doing a dense factorization of a 2D problem. – Nick Alger Nov 20 '16 at 20:36

1 Answers1

5

The "100,000 unknowns" rule-of-thumb applies to sparse matrices rather than dense ones. A naive direct solver, which doesn't take advantage of sparsity at all, could in principle have $O(n^3)$ complexity for some matrices. In practice, a good direct solver will have much lower computational complexity.

A better starting point would be to look at banded matrices. Factoring a matrix of bandwidth $d$ requires $O(d^2n)$ time, and the factors have bandwidth $d$ also. If you can reorder a sparse matrix so that it has small bandwidth, then you can make direct factorization much more manageable. Unfortunately, computing the optimal reordering is NP-complete. Fortunately, there are some good heuristics, i.e. the Cuthill-McKee algorithm, for reducing the bandwidth of a graph.

Now most FEM matrices, especially in 3D, actually have pretty high bandwidth. Rather than take this indirect approach of reducing fill-in by reducing bandwidth, you could ask which ordering, among the $n!$ possibilities, gives the least fill-in for a direct matrix factorization. This problem is also NP-complete, but there are also some good heuristics. These heuristics are the basis for approximate minimum degree orderings. A good reference is Amestoy's original paper on AMD ordering and Tim Davis's book.

As Bill Greene suggests, you should try it for yourself as well. scipy has a sparse LU factorization routine where you can specify to either leave the matrix in its natural ordering, or to use AMD. You can always try it out on the SuiteSparse matrix collection.

Daniel Shapero
  • 10,263
  • 1
  • 28
  • 59