Solve a large-scale linear system of equations with millions of unknowns

Question

I have a large-scale system of linear equations: $Ax = b$, where $A$ is an $n\times n$ square symmetric positive definite matrix (not sparse), $b$ is an $n \times 1$ vector and $x$ is $n\times 1$ unknowns. The magnitude of $n$ is several million.

1. How much memory does it take to store this problem, say, I use C++ and store numbers by double-precision?
2. How long does it take to solve it if I use a "reasonably good" computer and C++ language with a well-known library?
3. Any library (and its package/function) and algorithm to recommend for this problem?

Again, $A$ is an $n\times n$ square symmetric positive definite dense matrix.

(Further questions appended below) Thanks for the answers. I have further questions:

1. Matrix $A$ can be generated from several much smaller matrices whenever computation need its entries, then memory usage can be much smaller. Does this help to shorten the time and demand for hardwares?
2. If I only need to get the solution in 3 months, looks like using 1 GPU is enough. Am I right on this estimation?
3. The computation I do need to do very often is to take the solution of $x$ and multiply it with another vector of the same size, say 5 million by 1. How long does it take if I only use 1 GPU?
4. Any recommended books on using GPUs to solve linear equations?

Thank you very much!

Answer 1

You're going to need a large cluster or a supercomputer to solve this.

Memory usage is like arc_lupus commented, a double-precision float takes 8 bytes and there will be 1e-6^2 entries. We store just half because of the symmetry, but that's still 4 TB.

The algorithm you'll want is either a Cholesky factorization or Conjugate Gradient (CG). Choleksy is a standard routine in dense linear algebra packages. CG is usually used for sparse problems, so you might need to implement it yourself, but there's a possibility it'd be faster. As the comments allude to, you'll need distributed computing for both algorithms. You're problem may also benefit from techniques such as hierarchical low-rank structures. The libraries I know of include the following (They all have a Cholesky routine and the building blocks to implement CG):

• ScaLAPACK - the classic distributed LA library. CPU only. Intel MKL and IBM ESSL provide implementations with more optimizations.
• SLATE - CPU and GPU support
• DPLASMA - There are experiments with GPUs and hierarchical low-rank Cholesky, but I don't know if they've made it into the master branch yet.
• Elemental - Another CPU code. It appears LLNL is developing a GPU-accelerated fork, but I don't know the status of that.

My research involves dense non-symmetric problems, and I can solve a problem of size $n=256000$ on 48 NVIDIA V100 GPUs (16 GB of memory each) in about 75 seconds with the SLATE library; an SPD problem should be doable in about that time with half the GPUs. (I expect other GPU-accelerated libraries to perform roughly similar). Because memory scales as $O(n^2)$ and computation scales as $O(n^3)$, I expect you could solve the problem in about $75\times 4 = 300$ seconds on 384 GPUs.
I'm not sure on CPU solver performance, but you'll probably need thousands or tens-of-thousands of cores depending on how long you're willing the solve to take.