I endorse the recommendation by #stevengj to use preconditioned conjugate gradient. Computational efficiency then boils down to finding a good preconditioners. Given the finite element source of the problem, I expect multigrid to do well. Algebraic multigrid offers an easy gateway to multigrid. Not sure how versatile the implementations in GitHub - JuliaLinearAlgebra/AlgebraicMultigrid.jl: Algebraic Multigrid in Julia are. Alternatives are in e.g. GitHub - JuliaGPU/AMGX.jl (on GPU) and GitHub - JuliaParallel/PETSc.jl: Julia wrappers for the PETSc library (look for AMG under preconditioners). Other implementation surely exist.
Good luck. Please keep us posted on what you learn.