I’ve been facing a problem of solving linear systems Ax=b arising from discretized PDEs (Stokes equations in particular). Nively, it seems that solving Ax=b should not take much more time than simply multiplying matrix A by some vector v as long as A is sparse. More specifically, when discretizing a linear PDE, one end up with a sparce matrix A whose i:th row contains as many non-zero elements as there are nearest neighbors of node i within the spatial grid, and this is typically a small number. Thus, it seems that solving such a system should require O(n) time, where n is the size of the (square) matrix A. However, experience indicates that solving a linear system is way more costly than performing matrix multiplication, even in the case when A is sparse. In practice, I use MUMPS solver (that uses LU-factorization) implemented in PETSc library. I find that explicit solution becomes essentially unfeasible when the size of the matrix exceeds 100,000 or so. On the other hand, performing simple arithmetic manipulations with 100,000 floats seems like a trivial task. I come to the same conclusion experimenting with sparse matrices in MATLAB. For example:
A = sprand(100000,100000,0.0001) + speye(100000,100000);
v = ones(100000,1);
Now, multiplying A by v takes 0.02 seconds. However, solving A\v takes a much larger amount of time, even though A has roughly 10 elements on each row.
Is there a simple reason why solving a sparse system appears much harder than performing matrix-vector multiplication? Is there some solution to this?