Hey there,
I'm using Eigen to simulate fluids like water or smoke in a grid-based (eulerian) setup.
Currently I'm trying to simualte a 3D grid with 128*128*128 cells but get an allocation error when filling the vector or setting the matrix from triplets.

I am creating a Sparse Matrix of size (128^3, 128^3) to hold interpolation weights which I then use to evaluate the values in grid cells for the next timestep.

The basic idea is to follow an imaginary fluid particle from the cell center of a cell with 1D index i back in time to the position where it came from and to get the linear interpolation weights of the surrounding cells (max 8 at this time) and add them to the corresponding row with index i.
Unfortunately in a second step I forward trace an imaginary particle for cells that are in a specific condition, to add weights to the rows of other cells.
Therefore I am not able to tell how many non-zeros per row/column I will have in advance.

My current approach is to collect all interpolation weights in a triplet list (std::vector<Triplet>) and to then fill the matrix from that vector.
This works well for medium sized grids up to 100^3 but fails for a grid of size 128^3 in later simulation steps, where a lot of interpolations weights are included in the sparse matrix.

The program occupies around 500 MB for simulation relevant data in the 128^3 case.
When I fill the vector with Triplets it becomes quite huge and occupies an additional ~400 MB. For some reason, even though the program doesn't even exceed the 1GB (with 2GB maximum since I am running 32bit) I either get std::bad_alloc when the vector tries to allocate more memory, or if that does not happen, when I create the sparse matrix from triplets.

As you can imagine, I'm also concerned with the speed, so just inserting (insert/coeffRef) Elements directly into the matrix without being able to approximate the amout of non-zeros per row/column is not really an option.

So my question is: Is there any memory limitation on the sparse matrix (does it try to allocate a contiguous region of memory) - if so, how can I get around that? I assume that the vector memory allocation is the same problem, since it guarantees to reserve a contiguous region.
For me it is not important whether data is stored as whole or as fragments, I just need to be able to use the whole 2GB I have.
Any help is appreciated.
Thank you!

Relevant snippets:

Code: Select all

// create triplet list
std::vector<Triplet> tripletList;


then fill it while looping over all cells, for example

Code: Select all

tripletList.push_back(Triplet(currCell, uniqueIndizes[l], weightValues[l]));


where currCell and uniqueIndizes[l] are both 1D-Indizes for a grid cell -
and finally fill the matrix:

Code: Select all

// I use typedef Eigen::SparseMatrix<double> SpMatColMajor;
SpMatColMajor A(in->getCellCount(), in->getCellCount());
A.setFromTriplets(tripletList.begin(), tripletList.end());


Hi,

> get an allocation error when filling the vector or setting the matrix from triplets
Is it really during filling the matrix or during filling the vector?
How did you construct the sparse matrix? Are the data types capable for storing the values ?

> This works well for medium sized grids up to 100^3 but fails for a grid of size 128^3 in later simulation steps
This could depend on the requirement that vectors need contiguous memory.
I had a similar problem with very big vectors, a workaround could be to use "single linked lists" (e.g. std::forward_list) instead of a vector.
Don't use std::list it consumes to much memory!

> Is there any memory limitation on the sparse matrix (does it try to allocate a contiguous region of memory) - if so, how can I get around that?
Try for big matrices:
typedef SparseMatrix<Scalar, ColMajor, long long> SpMat;
>I assume that the vector memory allocation is the same problem, since it guarantees to reserve a contiguous region.
See work-around suggestion!

Best
Martin

Thank you for your quick response!

M00nMan wrote:Hi,
> get an allocation error when filling the vector or setting the matrix from triplets
Is it really during filling the matrix or during filling the vector?
How did you construct the sparse matrix? Are the data types capable for storing the values ?

Up to now I use

Code: Select all

typedef Eigen::SparseMatrix<double> SpMatColMajor;

M00nMan wrote:> This works well for medium sized grids up to 100^3 but fails for a grid of size 128^3 in later simulation steps
This could depend on the requirement that vectors need contiguous memory.
I had a similar problem with very big vectors, a workaround could be to use "single linked lists" (e.g. std::forward_list) instead of a vector.
Don't use std::list it consumes to much memory!

I now use

Code: Select all

std::forward_list<Triplet> tripletList;
// and for filling...
tripletList.emplace_front(Triplet(currCell, uniqueIndizes[l], weightValues[l]));


That seems to solve the problem for the triplet list (even though deallocation of the std::forward_list at the end of the method after matrix filling seems to be terribly slow), and I am now able to store more triplets (it did not fail to allocate now)!
But the std::bad_alloc still appears, now when I fill the matrix from the triplet list:

Code: Select all

SpMatColMajor A(in->getCellCount(), in->getCellCount()); // where getCellCount is 128*128*128
A.setFromTriplets(tripletList.begin(), tripletList.end()); // where tripletList ist a std::forward_list


So I still assume that the Sparse Matrix tries to reserve a contiguous memory block that is too big (since again I have enough memory so that it's not likely that I really am out of memory).
Any more suggestions on how to solve this problem?

M00nMan wrote:> Is there any memory limitation on the sparse matrix (does it try to allocate a contiguous region of memory) - if so, how can I get around that?
Try for big matrices:
typedef SparseMatrix<Scalar, ColMajor, long long> SpMat;
>I assume that the vector memory allocation is the same problem, since it guarantees to reserve a contiguous region.
See work-around suggestion!

Best
Martin

This unfortunately doesn't work since in PlainObjectBase.h I get a warning (which is considered an error in my project) in the check_rows_cols_for_overflow while trying to evaluate the max_index : warning C4293: '<<' : shift count negative or too big, undefined behavior eigen\src\Core\PlainObjectBase.h 42

But the indexing with int should not be a problem, so I don't assume that this is necessary anyway.

Thank you!
Best,
Niklas

>That seems to solve the problem for the triplet list (even though deallocation of the std::forward_list at the end of the method after matrix filling seems to be terribly slow), and I am now able to store more triplets (it did not fail to allocate now)!

Deallocate the first element on the fly direct after an insert! This has twice advantages first you'll have more free space second deallocation is faster.

>So I still assume that the Sparse Matrix tries to reserve a contiguous memory block that is too big (since again I have enough memory so that it's not likely that I really am out of memory).
>Any more suggestions on how to solve this problem?
Iterate over the std::forward_list and count the NNZ per column than allocate the sparse matrix. This is much faster because reallocation and thus resorting during insert it's not necessary. Maybe the could also avoid a std::bad_alloc.
A.reserve(NNZs);

>But the indexing with int should not be a problem, so I don't assume that this is necessary anyway.
This depends on the number of NNZs as far as I know the count is also represented with the Index data-type.

>This unfortunately doesn't work since in PlainObjectBase.h I get a warning (which is considered an error in my project) in the check_rows_cols_for_overflow while trying to evaluate the max_index : warning C4293: '<<' : shift count negative or too big, undefined behavior eigen\src\Core\PlainObjectBase.h 42
I don't think this warning depends on the typedef.

I use the following in my project:

Code: Select all

class spmatrix : public Eigen::SparseMatrix<double,Eigen::RowMajor, long long int> {
public:
        typedef Eigen::SparseMatrix<double,Eigen::RowMajor, long long int> Base;
...


This works fine!

--
Best
Martin

M00nMan wrote:Deallocate the first element on the fly direct after an insert! This has twice advantages first you'll have more free space second deallocation is faster.

To be honest I don't really understand this point. Could you give a minimum example on that please?

M00nMan wrote:Iterate over the std::forward_list and count the NNZ per column than allocate the sparse matrix. This is much faster because reallocation and thus resorting during insert it's not necessary. Maybe the could also avoid a std::bad_alloc.
A.reserve(NNZs);

That seems to be exactly the problem. This is what Eigen does when you call setTriplets anyway - they count the number of non zeros per column (since I'm using column major matrix). I just ran the example wherre i get bad_alloc and counted the number non zeros beforehand. The max number is 28 so the Sparse matrix will try to allocate 28 * 128 * 128 * 128 * 16 =~ 940 MB for the matrix (even though most columns will have far less than that) so the bad_alloc has a point since that would get me somewhere around 2GB memory.
However the list only contains 16761078 elements, so it has a (pure storage size) of 16761078 * 16 =~ 270 MB.

So the problem really is that the Sparse matrix will be huge in the uncompressed mode..
This probably leaves me with simulating smaller grids or directly inserting into the matrix I guess?

>To be honest I don't really understand this point. Could you give a minimum example on that please?

In set_from triplest ther is something like:

Code: Select all

    for(InputIterator it(begin); it!=end; ++it)
      trMat.insertBackUncompressed(it->row(),it->col()) = it->value();

After some function modification:

Code: Select all

    for(InputIterator it(begin); it!=end; ++it){
      trMat.insertBackUncompressed(it->row(),it->col()) = it->value();
  // protect the iterator !
    if(it != begin)
      mylist.pop_front();
      }

>The max number is 28 so the Sparse matrix will try to allocate 28 * 128 * 128 * 128 * 16 =~ 940 MB for the matrix (even though most columns will have far less than that) so the bad_alloc has a point since that would get me somewhere around 2GB memory.
>However the list only contains 16761078 elements, so it has a (pure storage size) of 16761078 * 16 =~ 270 MB.

Count the NNZ per row / col and use instead the max value of NNZ a vector of NNZs for each row / col like:
A.reserve(NNZs_vect);
Than exactly 16761078 NNZ plus some over head for indexing will be allocated. I assume this should not more than 350 MB (without guarantee ).
You can take a look in in set_from_triplets there this is done.

>So the problem really is that the Sparse matrix will be huge in the uncompressed mode..
That a point of definition; my matrices are 450 GB ~ 670GB with 22.000.000.000 elements of type double.

--
M00nMan

I get the bad_alloc when I call reserve on the matrix with the vector of NNZs (basically it doesn't really matter if I do it beforehand oder if it happens when set_from_triplets is called)
therefore I can't even use your approach to remove the element from the forward list after inserting it into the matrix.

I don't really understand why this is not working if the Sparse Matrix really is only using the real number NNZs per column.

>I don't really understand why this is not working if the Sparse Matrix really is only using the real number NNZs per column.
Try the sparse matrix typedef with data type long or long long.

>I get the bad_alloc when I call reserve on the matrix with the vector of NNZs
A sparse matrix consists of vectors may be there is not enough contiguous memory available.

BTW: Do you use the Eigen parallel option?
I had some trouble during big matrix allocation when it was used.