Since there are many people in this forum that have a wide knowledge in not only programming but also in numerical analysis, I would like to ask for a recommendation on which solvers to consider to solve stardard linear systems Ax = b, where A is a very large sparse matrix (around 3GB of RAM). This matrix has the property of being positive definite and symmetric.
I have seen that LinearSolve.jl has many methods, but I feel the documentation does not give advice towards which solvers should be considered in the case that I am interested.
I must also mention that I do not want to solve the same linear system for different b, since the A matrix will be changing from one solve to another.
Try sparse Cholesky cholesky(A) \ b first (the built-in one from CHOLPACK in SparseArrays.jl, or you could try Pardiso.jl or others). (For a sufficiently big problem, this may run out of memory, depending on the sparsity pattern.) Sparse-direct solvers like this have the advantage of being relatively bulletproof with no adjustable parameters: they either work, or they run out of memory.
Also try conjugate gradient (this is usually the best iterative method for SPD matrices). May require some tuning of the tolerance/convergence parameters, depending on the matrix and how accurate you need the solution to be.
Ideally use preconditioned conjugate gradient, but that requires you to figure out a preconditioner, which is problem-dependent and requires some knowledge and experimentation. RandomizedPreconditioners.jl is one option that is easy to try.
Thank you very much for the suggestions. I will comment my findings.
As suggested by @stevengj, I first tried doing sparse Cholesky, but my real problem ran out of memory (I tried this on a smaller problem and it worked beautifully). The conjugate gradient methods were taking too much time, probably because I was not using a suitable preconditioning.
As per your suggestion, I tried the Algebraic Multigrid methods (which I did not know anything about, so thank you for pointing them out) and I had incredibly good results. Using the example here from AlgebraicMultigrid I found the solution of my system extremely quickly (less than 5 seconds for my system, which as I have said is quite large and occupies 3GB of memory in sparse form).
So I will spend a little bit learning about these methods but they are clearly what I need, so I do not know if you have any additional suggestion but thank you again for pointing me in this direction.
I mean the space dimension of the domain discretized with the finite element method. If it is 2 , sparse direct solvers and cholesky factorization are still competitive. If it is 3, the fill-in (reduction of sparsity during the factorization) from these methods scales badly.
What your revealed so far speaks for space dimension 3.
You als may test AMGCLWrap.jl for a preconditioner which wraps an algebraic multigrid library written in C++.
