There are many different linear solvers, some which work best for diagonally dominant matrices, some for symmetric, some for positive definite ones, some for banded matrices, etc... There are direct methods, iterative methods, Krylov space methods, etc... Is there a good heuristic to choose a linear solver that scales well in parallel as the problem size increases? In particular, is there some sort of block diagram that can specify an efficient (or optimal) solver for specific matrix system types?
$\begingroup$
$\endgroup$
4
-
1$\begingroup$ Similar question: scicomp.stackexchange.com/questions/81/… $\endgroup$– David KetchesonJan 11, 2012 at 19:15
-
$\begingroup$ David: I think the difference is "sparse" versus general linear solvers. $\endgroup$– aeismailJan 11, 2012 at 21:39
-
$\begingroup$ I don't see the difference myself, unless people are using Krylov for dense systems now and I'm just out of touch... $\endgroup$– J. M.Jan 11, 2012 at 23:55
-
$\begingroup$ @DavidKetcheson: I agree, it's broad, but I think it's fair to ask for guidelines on first steps towards choosing a linear solver before delving into the weeds and looking at what is the absolute best algorithm to choose for a given problem. $\endgroup$– Geoff OxberryJan 12, 2012 at 19:58
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
I presume that by "optimal" you mean "fastest solve".
I think the PDE solver guys have the best advice when it comes to advice about the most appropriate algorithm for certain discretizations and problems, such as Helmholtz problems (see, for instance, this question) and other cases where problem structure can be exploited.
From a basic numerical linear algebra perspective, here are some general guidelines:
If your problem is dense, not low-rank, and small enough to fit into memory, direct solvers are a decent choice.
If your problem is ill-conditioned, preconditioning or regularizing is a good idea.
If you can exploit structure (bandedness, symmetry, diagonal dominance), doing so will save you operations.
If your problem is sparse and small enough, sparse direct solvers can work well.
If your problem is sparse and large, use an iterative solver, because fill-in will render sparse direct solvers useless. The efficacy of an iterative solver is problem-dependent (see this question) and depends on spectrum. You will probably want to try all of them, unless your problem has a specific structure, or you know from experience that a certain algorithm works well for your problem.
If your problem is dense, and low-rank, you may be able to exploit that low-rank structure, either through use of iterative solvers, or specialized algorithms (which require special problem structure, as alluded to above with PDEs).
When it comes to parallelism, the best algorithms minimize communication, which is accomplished by structuring both the algorithm and the data accordingly. I am not an expert in these matters. I know that block cyclic matrix decompositions are popular for dense numerical linear algebra, and...that's it. I work almost entirely on serial algorithms, so I leave the discussion of parallelism to others.
answered Jan 12, 2012 at 9:13
$\begingroup$
$\endgroup$
1
It's probably not even close to optimal, but here's the basic hueristic I use:
Don't use a direct solver unless the matrix involved is diagonal, tridiagonal, can be easily put into tridiagonal form (as with the split operator method), or has some other property that makes it easy to directly solve.
Use whatever iterative solver your favorite library uses otherwise. Iterative solvers tend to be nicer to parallelize.
Note that there are lots of specific cases where you really need some specific algorithm, as indicated by this answer.
-
3$\begingroup$ Your answer seems to be overlooking sparse direct methods, which are essentially $O(n^{3/2})$ complexity in 2d and $O(n^2)$ complexity in 3d for discretizations with local basis functions. $\endgroup$ Jan 13, 2012 at 14:43