Direct solvers are your ‘black box’ for most of your needs.
In infinite precision, direct solvers would find the exact solution.
The algorithm has a fixed number of steps, therefore no tolerances need to be used.
They are generally the most robust algorithms for ill-conditioned systems.
They scale terribly (both in system size and parallelization). Unless you have a specific performance problem, you should default to a direct solver, but be aware of the limits.
Sparse matrices are your saving grace! Do your best to protect them (which is why you should always store the factors, and NEVER directly calculate the inverse unless you strictly need the elements of the inverse itself).
Modern LU methods can benefit dramatically from distributed and GPU computing, however this doesn’t completely eliminate the scaling issue.
Generally speaking, you won’t see a significant speedup with parallelization until you get to a very large number of compute nodes.
If you use them on a massive problem, start with a small sub-system and work your way upwards to test scaling.
Warning: Some standard implementations (like
numpy.linalg.solve) are sophisticated enough to handle singular matrices as well as non-singular ones under the hood. Be extremely careful and always verify your answers!