This question has also been asked at Stack Overflow.

This is a very general question regarding the maximum size of a set of linear equations to be solved by today's fastest hardware, in the form:

$$X = AX + B$$

Where we solve for $X$, and

  • $A$: $N\times N$ a sparse matrix of floats;

  • $B$: $N$-vector of floats.

This becomes $X(I-A) = B$ which is best solved using factorization (and not matrix inversion) as I read here.

Do you know or have a reference to a benchmark or paper which gives some maximum value for $N$ with today's fastest hardware? Assume that I can have that hardware, this is a sort of mental what-if exercise. Most benchmarks I have seen use $N < 10000$. I am thinking about $N>10^7$ or more to be processed within a month.

Please consider not only the computational dimension but also storage for $A$. It can be a problem, e.g., assuming $N = 10^6$, storage would be $4\times 10^{12}$ bytes $\approx$ 4 Terabytes for a totally dense matrix, which is about manageable I guess.

Lastly, can the method to solve the system be parallelized so that I can make the assumption that with parallelization $N$ can become pretty large?

later addition to question:

One possible application of this is: consider N products from factories (e.g. screws, an electric motor, an engine, a car). Most products require parts from other products and are used as parts for other more complex products in addition to some quantities ending in the shops - e.g. an electric motor.

For the $i^{th}$ product we have: $X_i = a_{i,1} X_1 + a_{i,2} X_2 + \cdots + a_{i,N} X_N + B_i$ meaning that the quantity of $X_i$ we aim to produce requires $a_{i,1} X_1$ items from product $X_1$, to $N$, plus $B_i$ which is what ends up in shops for general consumption.

Our aim is to find the solution to this problem given final demands ($B$) and "inter-dependancy" quantities $A$.

I guess $A$ is sparse but this probably depends on applications, e.g. Toyota states that a basic car consists of 30,000 other products (all broken down to screws). This is why I mention sparse but I also go for the worst case of "mild-denseness".

  • 2
    $\begingroup$ Do you have a specific scientific problem in mind that leads to such a linear system? (Idle speculation might be interesting to you, but difficult to answer for us...) For practical purposes, you should look at Krylov methods; if you want to know what today's fastest hardware is capable of, you can look at the results of the current LINPACK benchmarks, which involves solving a (dense) linear system. $\endgroup$ Sep 25, 2015 at 17:06
  • $\begingroup$ i have added an example application in the question. thanks $\endgroup$
    – bliako
    Sep 26, 2015 at 11:52

1 Answer 1


For your applications, dense matrix techniques such as the different factorization techniques will not work. Plus your matrix is sparse anyway so you should be using an iterative solver (e.g. Krylov subspace methods). These easily handle problems of the scale you are interested in and there are several sparse matrix solver libraries implemented for parallel architectures (e.g. Petsc, Trilinos, others).

There are several advantages to using sparse matrix techniques:

  1. Storage. Who wants to store a dense matrix? Do you have terabytes of RAM waiting to be used? Most people don't, and you should not want to have to use disk space because that will cripple any sort of algorithm performance. You've already stated that your matrix A is sparse, which generally means that total storage space should only be $O(N)$. So it should only take megabytes to store your matrix and required vectors, not terabytes.

  2. Time. Dense matrix factorizations take $O(N^3)$ operations, even for sparse matrices. On a serial processor, a system of equations with $N>10^7$ can easily take on the order of years to compute a factorization. Even when parallelized you still have to deal with storage and communication costs. Sparse methods rely on matrix-vector multiplication, which is an $O(N^2)$ operation for dense matrices but a $O(N)$ operation for sparse matrices. And usually it requires less than $O(N)$ iterations to converge. This will take much less time in general (and can take even less time when a good preconditioner is used).

  3. Parallelization. Because iterative methods rely on matrix-vector multiplication and vector-vector addition, they are highly parallelizable (especially on shared memory architectures, but also distributed memory too). Many libraries already have parallel support (again see here) for various architectures.
  • 1
    $\begingroup$ A slightly nitpicking point: Iterative sparse methods depend on matrix-vector multiplications, but there are also direct sparse methods which do not and instead factor the matrix whilst taking advantage of the sparsity. These will also typically scale much less than O(N**3) for sufficiently sparse matrices. $\endgroup$
    – Ian Bush
    Sep 26, 2015 at 7:33
  • $\begingroup$ Thanks for your answer. It clarified things. I have now started using PETSc in order to experiment with size and denseness of the matrix A. If there are any pointers to tutorials or code (other than the standard ones on the PETSc website) or any book, please let me know. So far I have managed to write a program to produce random A and B with variable density, however applying a KSP solver proves much harder. Many thanks again. $\endgroup$
    – bliako
    Sep 28, 2015 at 11:32

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.