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!
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.