5

Problem: A is square, full rank, sparse and banded. It has way too many elements to be stored as a single matrix in Matlab (at least ~4.6*1018 and ideally ~1040, both of which exceed max array size. EDIT: A is stored as sparse, and the problem is not with limited memory but with limited number of elements). Therefore I have to store it as a collection of smaller arrays (rows/diagonals/columns/blocks).

Looking for: a way to solve Ax=b, with A given as a collection of smaller arrays. Ideally in Matlab but not a must.
Alternatively, if not in Matlab: maybe there's a program that can store and solve such a big A?

Found so far: methods if A is tri/pentadiagonal, but my A has N diagonals. Also found something about partitioning A to blocks, but couldn't find a way to then solve a linear system with these blocks.

p.s. The system is 64-bit.

Thanks everyone!

OOO
  • 51
  • 2
  • 2
    Have defined your matrix as a normal matrix? Not explicitly as [**sparse**](http://www.mathworks.de/de/help/matlab/ref/sparse.html)? The latter should save you a lot of memory. For solving: [**this article**](http://www.mathworks.de/de/help/matlab/math/sparse-matrix-operations.html) is probably helpful. I've never used sparse matrices and can't help you further, but surely others can ;) – Robert Seifert Sep 14 '14 at 17:20
  • `sparse` has size limits as well. Check the second output argument of `computer`, it returns the maximum number of elements that can be indexed. `[~,maxsize]=computer`. I don't link them here because they depend on your matlab version, but as far as I know no Version supports 4.6*10^18 elements. – Daniel Sep 14 '14 at 18:08
  • You need to find a solver that will use a function pointer/handle rather than an array stored in memory. I just quickly looked but didn't find anything. However, I am sure one exists. With that being said, I suspect you will have precision issues with that many elements regardless of the solver. – AnonSubmitter85 Sep 14 '14 at 20:53
  • AnonSubmitter85: Sounds promising, do you have any keywords or search directions for this kind of solver? – OOO Sep 14 '14 at 21:35
  • @OOO. Sorry for the comment earlier, i didn't catch that the limiting factor was the number of elements and not the memory. As Anon suggested, for this size of system you will have to find some special algorithms. As a starter you can have a look at this method [Fast Inverse using Nested Dissection](http://mc.stanford.edu/cgi-bin/images/0/04/Li_phd.pdf). It's a long read but promising. I don't know if there are ready made solver using this algorithm yet. – Hoki Sep 15 '14 at 08:58
  • @OOO, what number of elements are you talking about? The number of nonzero or the the total number of elements of the full matrix? What is the number of row of the associated full matrix? When you said N diagonals, what is the value of N? – innoSPG Jul 09 '15 at 18:52

3 Answers3

1

Not using Matlab would allow you to store larger arrays. ROOT is an open source framework developed at CERN that has C++ and Python interfaces and a variety of solvers. It is also capable of handling huge datasets and has a variety of visualization and analysis tools as well.

If you are interested in writing C or Fortran BLAS(Basic Linear Algebra Subroutines) and CBLAS would be good options. There are many open source and proprietary implementations of BLAS that should be available for most Linux/UNIX distributions. There are also plenty of examples showing how to use the BLAS subroutines in C and Fortran code available online.

Matt
  • 545
  • 3
  • 16
0

If you have access to MATLAB's Parallel Computing Toolbox together with MATLAB Distributed Computing Server, you may be able to store A as a distributed array, in other words a single array whose elements are distributed across the memories of multiple machines in a cluster. You can call MATLAB's backslash command directly on a distributed array, and MATLAB handles the parallelization for you.

Sam Roberts
  • 23,951
  • 1
  • 40
  • 64
  • Sounds interesting thanks! I looked them up and couldn’t find the following: Do you know if distributed arrays can contain more elements than regular arrays? I thought the problem is that the elements index is limited. Do distributed arrays overcome this? – OOO Sep 14 '14 at 21:33
  • @OOO The distributed array may be able to contain more elements than a regular array **if** your different workers are on **separate** machines (so you actually increase your total memory available). Otherwise distributing work on several cores of the same machine will only gain you processing time, not extra memory. – Hoki Sep 14 '14 at 23:54
  • @Hoki Thanks. Just to make sure: my problem is not with memory but with number of elements (i.e., an empty sparse matrix has a limited size as per this page: http://www.mathworks.com/matlabcentral/answers/91711-what-is-the-maximum-matrix-size-for-each-platform), will distributed array on multiple machines have an increased number of elements? – OOO Sep 15 '14 at 00:47
0

I wanted to put this as a comment, but I think it is better to state it as an answer. You have a serious problem. It is not only a problem of indexing, it is also a problem of memory: 4.6x10^18 is huge. That is 4.6 exa elements. If you store them as real single precision, you need 4x4.6 exabyte of memory. A computer which such a huge memory, does not yet exists to my knowledge. You will need to gather all the storage (hard disk, not RAM) of a significant proportion of all computers in the world to store such a matrix. Think about it. Going to 10^40 elements is nearly impractical for the time being. With your 64 bit computers, the 64 bit address space can bearly address 4.6x10^18 elements. 64 bits address (or integer) makes it possible to directly index 2^64 elements which is roughly 16x10^18. So you have to think twice.

Going back to the problem itself, there are chances that you can turn your matrix into an implicit operator. By implicit operator, I mean, you do not need to store it, because it has a pattern that you know how to reproduce, or you can apply it to a vector without actually forming the matrix. If you have the matrix in hand, you are very likely in this situation, considering what I said above. If that is the case, to solve your problem, you simply need to use an iterative solver and provide a black box that does your matrix multiplication. Going to other directions might be a waste of your time.

innoSPG
  • 4,588
  • 1
  • 29
  • 42
  • OP has stated that it's a sparse matrix. – rlbond Jul 09 '15 at 00:52
  • @rlbond I do not get your point. Can you explain a little bit? – innoSPG Jul 09 '15 at 00:54
  • @rlbond, I still do not get it. I have some experience in sparse matrices. Can you point to something I said that is in contradiction with what OOO said, or in contradiction with sparse matrices? – innoSPG Jul 09 '15 at 01:09
  • OP said their matrix is mostly zeros, those do not use memory in a sparse matrix data structure. – rlbond Jul 09 '15 at 17:55
  • I got you now. OOO also said that their matrix is already stored as sparse. From that I deduced that he was talking about the number of nonzero. I'am putting a comment to ask him for clarification. – innoSPG Jul 09 '15 at 18:50