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!

  • 6
    $\begingroup$ 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. $\endgroup$ Nov 20 '16 at 10:26
  • 4
    $\begingroup$ 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. $\endgroup$ Nov 20 '16 at 11:08
  • 2
    $\begingroup$ @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: scicomp.stackexchange.com/questions/3262/… $\endgroup$ Nov 20 '16 at 14:43
  • 2
    $\begingroup$ 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. $\endgroup$ Nov 20 '16 at 14:54
  • 1
    $\begingroup$ 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. $\endgroup$
    – Nick Alger
    Nov 20 '16 at 20:36

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.


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