# Solving huge dense square symmetric linear system

I have a linear system of the type

$$A x = y$$

where A is a dense, square, symmetric, positive definite matrix, $$x$$ a vector of unknown parameters, and $$y$$ is a vector of observed quantity. I know that using a Cholesky decomposition based algorithm should be the fastest way to find the values of $$x$$ that fulfil the equation. The only problem is I would like to do so for a matrix $$A$$ with dimension 500000x500000, and $$x$$ and $$y$$ 500000x1. Is there any scientific computing library to do so? Which kind of hardware should I use?

• Well I can tell you it is doable - on 16384 cores of a modern cluster a collaborator and I have demonstrated finding all evals and evecs of a dense, real, symmetric matrix of about that size which is a somewhat tougher calculation than you need. We used ScaLAPACK which also has solvers in, being the distributed memory version of LAPACK. So it is possible, but whether ScaLAPACK is a good way to go for your problem I feel less sure - hence a comment rather than an answer (basically I do diagonlisation, I rarely solve equations) Jan 19, 2022 at 14:05
• You're going to need a system with a lot of memory to even store the matrix A. Multiply 500000 by 500000 by 8 bytes per entry and you'll have 2000 gigabytes of storage required for A. This isn't something you can do on a desktop machine, but could well be within the capacity of a high performance computing cluster. Jan 19, 2022 at 15:28
• But even if possible, it seems questionable whether whatever you are doing is really best described by a dense matrix! Jan 19, 2022 at 21:20
• @WolfgangBangerth why questionable? I can think of numerous problems that are really best described by a dense matrix. Jan 20, 2022 at 9:10
• @BrianBorchers You're missing a division by 2 because $A$ is symmetric (but that doesn't change the conclusion). Jan 20, 2022 at 12:34