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!

  • 2
    $\begingroup$ You can estimate the memory size by assuming double = 8 bytes. The vector size is negligible compared to the matrix. When assuming n = 1e6 you'll look at a matrix with a size of 1e12 entries, corresponding to (at least) 8 TByte of memory. Not sure if it is possible at all to solve such a system within reasonable time on a "reasonably good" computer. $\endgroup$
    – arc_lupus
    Commented Feb 21, 2022 at 8:03
  • $\begingroup$ It's doable - see my comment to scicomp.stackexchange.com/questions/40880/… . You just need access to a cluster with at least a few thousand cores. $\endgroup$
    – Ian Bush
    Commented Feb 21, 2022 at 9:30
  • $\begingroup$ Some more details about the underlying problem that generates this large dense matrix might help someone provide a much more efficient means of solving this system. $\endgroup$ Commented Feb 22, 2022 at 14:55
  • $\begingroup$ you can also use "row-action" methods which iterative solve the system for small subsets of rows, this can be orders of magnitude faster than working with whole matrix at each step $\endgroup$ Commented Sep 1, 2022 at 22:45

1 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.

  • 1
    $\begingroup$ Just curious, is that a typo or real: 48 NVIDIA V100 GPUs (16 TB of memory each)? Or should it be 16 GB? $\endgroup$
    – arc_lupus
    Commented Feb 21, 2022 at 14:04
  • 1
    $\begingroup$ That was a typo. It should be 16 GB. $\endgroup$ Commented Feb 21, 2022 at 15:32
  • $\begingroup$ Everything else would have been really interesting :-D. $\endgroup$
    – arc_lupus
    Commented Feb 21, 2022 at 16:04
  • $\begingroup$ Thanks for this interesting answer! I had always assumed that the largest scale linear systems were solved via some kinda krylov subspace method rather than just running a decomposition algorithm in parallel, but I guess methods like that, including CG, are best left to problems with specific structure which accelerates matvec computation? $\endgroup$ Commented Feb 22, 2022 at 19:48
  • $\begingroup$ Depends on the problem. CG can be useful if it converges fast enough to the desired precision. CG is bound by memory, and network performance w/ lots of synchronizations while Cholesky can get close to the peak arithmetic rate, so high iteration counts will make CG costly. Preconditioning is very helpful but requires tuning/knowledge of the problem. Cholesky is also useful for cases with multiple right hand sides. It's probably worth trying both if you going to be solving more than a couple systems, but people tend to look to Cholesky for dense problems. $\endgroup$ Commented Feb 23, 2022 at 1:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.