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?