12
$\begingroup$

I need to solve the same sparse linear system (300x300 to 1000x1000) with many right hand sides (300 to 1000). In addition to this first problem, I would also like to solve different systems, but with the same non-zero elements (just different values), that is many sparse systems with constant sparsity pattern. My matrices are indefinite.

Performance of the factorization and initialization is not important, but performance of the solve stage is. Currently I'm considering PaStiX or Umfpack, and I will probably play around with Petsc (which supports both solver) Are there libraries capable of taking advantage of my specific needs (vectorization, multi-threading) or should I rely on general solvers, and maybe modify them slightly for my needs ?

What if the sparse matrix is larger, up to $10^6 \times 10^6$ ?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
10
$\begingroup$

Without taking sides the discussion about whether to use direct or iterative solvers, I just want to add two points:

  1. There exist Krylov methods for systems with multiple right-hand sides (called block Krylov methods). As an added bonus, these often have faster convergence than standard Krylov methods since the Krylov space is built from a larger collection of vectors. See Dianne P. O’Leary, The Block Conjugate Gradient Algorithm and Related Methods. Linear Algebra and its applications 29 (1980), pages 239-322. and Martin H. Gutknecht, Block Krylov space methods for linear systems with multiple right-hand sides: An introduction (2007).

  2. If you have different matrices with the same sparsity pattern, you can precompute a symbolic factorization for the first matrix, which can be reused in computing the numerical factorization for this and the subsequent matrices. (In UMFPACK, you can do this using umfpack di symbolic and passing the result to umfpack_di_numeric.)

Improve this answer
$\endgroup$
9
$\begingroup$

There is typically a trade-off between the amount of work you put into constructing a good preconditioner for an iterative solver and the work you save by using a good preconditioner when actually solving the linear systems. In your case, the case is pretty clear: put as much work as you can into constructing a good preconditioner because you have to solve so many linear systems. In fact, I think it is appropriate to invest the time to get the perfect preconditioner: an LU decomposition (using UMFPACK, for example, or the Pardiso solver that comes as part of Intel's MKL). Then simply apply this decomposition as many times as necessary. If you have $O(N)$ linear systems to solve, nothing can be expected to beat an exact decomposition.

Improve this answer
$\endgroup$
3
  • 5
    $\begingroup$ Your last statement is debatable. Consider an exact multifrontal factorization of a 3D FEM or FD discretization over a cube, which should require $O(N^2)$ work and $O(N^{4/3})$ memory usage. The exact solves therefore require $O(N^{4/3})$ flops per right-hand side, and so, for sufficiently large $N$, any iterative solver with a lower asymptotic complexity will be faster. $\endgroup$
    – Jack Poulson
    Feb 8, 2013 at 15:46
  • 3
    $\begingroup$ Maybe. But as a matter of practical consideration, sparse direct solvers are still darn fast given that the constant in front of even an $O(N)$ solver is pretty large whereas the constant in front of the $O(N^{4/3})$ is not. $\endgroup$
    – Wolfgang Bangerth
    Feb 8, 2013 at 23:25
  • 3
    $\begingroup$ The catch is that you run of memory and patience for the factorization at about the crossover point. For the 7-point Laplacian, multigrid needs about 50 flops/dof, leading to a flops crossover (versus back-solve) at around $10^5$ dofs. The back-solve uses a lot more memory, but kernels for many right hand sides are commonly available. Multigrid is not usually written for many right hand sides, thus sacrificing vectorization potential. I'd wager that you can write an MG algorithm for which time-per-RHS is less than a CHOLMOD (or any other package) back-solve for 3D Laplacian at some $n<300k$. $\endgroup$
    – Jed Brown
    Feb 9, 2013 at 15:39
3
$\begingroup$

You're not quite clear in your statement of the problem when you talk about "the same non-zero elements (just different values)" Are you saying that the matrix has a constant sparsity pattern but the actual values change? Or, are you saying that the matrix is in fact constant?

Assuming that the sparse matrix is constant and only the right hand side is changing, then you should be looking at methods that use direct factorization (of the form $PA=LU$) of the matrix, and then solve for each right hand side by forward/backward substitution. Once the factorization is complete, each solution will be extremely fast ($O(n^2)$ time for completely dense factors but for sparse factors this will be proportional to the number of nonzeros in the factors.)

For multiple right hand sides and systems of equations of this size, iterative methods are typically not worth while.

All of the packages you mentioned offer direct factorization methods (although PetSc is mostly known for its iterative solvers.) However, your systems are so tiny that it is unlikely that you could get substantial parallel speedups, particularly in a distributed memory environment.

I'd suggest using Umfpack for this job- PaStix and PetSc are overkill.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks for your answer. In order to clarify: I asked first for a single matrix with many right hand sides, and then, another problem is a collection of matrices with the same sparsity patern but the values change, each of them must be solved for many rhs. Subsidary question: what if the sparse matrix is now 10^5x10^5 to 10^6x10^6 ? $\endgroup$
    – nat chouf
    Feb 8, 2013 at 16:43
  • 2
    $\begingroup$ My rule of thumb is that a sparse direct solver (taking as example the discretization of a 2d PDE) is faster than even good iterative solvers if the size of the matrix is less than $10^5$. That may be a rough guess, but it may give you an idea. $\endgroup$
    – Wolfgang Bangerth
    Feb 8, 2013 at 23:27
  • $\begingroup$ Using an iterative method for your larger systems with only a single right hand side might make sense, particularly if you don't need very accurate solutions and particularly if you can find an effective preconditioner or your systems are already well conditioned. However, if your systems are badly conditioned, you need accurate solutions, and you can't find a good preconditioner, then you'll likely still be better off with direct factorization. $\endgroup$
    – Brian Borchers
    Feb 9, 2013 at 2:29
  • $\begingroup$ Another important consideration is the memory requirements. Once you get into very large sparse systems where $N$ is of order $10^6$, you can easily run out of memory to store the direct factorization. At the point you'll definitely be forced to switch to using an iterative method. $\endgroup$
    – Brian Borchers
    Feb 9, 2013 at 2:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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