c++ - Very slow filling of sparse matrix and no memory gain in Eigen

Question

I am trying to implement an example that is given in the Eigen tutorial as pseudocode. As far as I understand, it illustrates the recommended method for filling a sparse matrix, provided that the number of non zero entries per column is known.

The pseudocode is found under the Header "Filling a sparse Matrix"and is written as follows:

1: SparseMatrix<double> mat(rows,cols);         // default is column major
2: mat.reserve(VectorXi::Constant(cols,6));
3: for each i,j such that v_ij != 0
4:   mat.insert(i,j) = v_ij;                    // alternative: mat.coeffRef(i,j) += v_ij;
5: mat.makeCompressed();                        // optional

My attempt to transform it into C is found below. I have (hopefully) written vee() such, that it will create 2500 non-zero elements per column. Therefore 2500 should correspond to 6 in the example. I set it to 3000 to test make.Compressed as well.

Unfortunately, I don't understand the behaviour of the programme. It does i=0...3000 in seconds, then is stuck for minutes. Then it goes to 6000 and is stuck again for minutes. Why is that and how to get better performance?

Furthermore, the memory usage is very strange. You can see that at times toward the end Eigen uses significantly more memory than needed for the corresponding dense matrix in GSL.The memory used also fluctuates wildly. in steps bigger than 100MB

I compile and run like this:

ludi@ludi-M17xR4:~/Desktop/tests$ g++ -o eigenfill.x eigenfill.cc -L/usr/local/lib -lgsl -lgslcblas && ./eigenfill.x

#include <iostream>
#include <stdio.h>
#include <stdlib.h>     
#include <Eigen/Sparse>
#include <gsl/gsl_matrix.h>
#define rows 1e4
#define cols 1e4



/*-- DECLARATIONES --*/
int FillMatrix(Eigen::SparseMatrix<double> mat);
double vee(int i, int j);

int main()
{
printf("---> Watch gsl matrix memory usage!
");    
gsl_matrix *testmat = gsl_matrix_calloc(rows, cols);
sleep(20);
gsl_matrix_free(testmat);
printf("---> Watch eigen matrix memory usage!
");  

Eigen::SparseMatrix<double> mat(rows,cols);         // default is column major
FillMatrix(mat);
printf("------------------------DONE");
return(0);
}

/*-- --*/
int FillMatrix(Eigen::SparseMatrix<double> mat)
{
int i, j;
Eigen::VectorXd Vres;
mat.reserve(Eigen::VectorXi::Constant(cols,3000));
    for(i=0;i<rows;i++)
    {
        if(i%500==0){printf("i= %i
", i);}
    for(j=0;j<cols;j++)
        {
        if (vee(i,j) != 0){mat.insert(i,j) = vee(i,j);    /*alternative: mat.coeffRef(i,j) += v_ij;*/ }
        }

    }
printf("--->starting compression"); 
mat.makeCompressed();

return(0);
}
/*-- --*/
double vee(int i, int j)
{
    double result = 0.0;

    if(j%4 == 0){result =1.0;}

    return result;
}
/*-- --*/

EDIT

One answer reminded me I needed to use addresses, because the local variables of FillMatrix() will be gone when it runs through. I tried the following, which does not compile:

#include <iostream>
#include <stdio.h>
#include <stdlib.h>     
#include <Eigen/Sparse>
#include <gsl/gsl_matrix.h>
#define rows 1e4
#define cols 1e4



/*-- DECLARATIONES --*/
int FillMatrix(Eigen::SparseMatrix<double> & mat);
double vee(int i, int j);

int main()
{
printf("---> Watch gsl matrix memory usage!
");    
gsl_matrix *testmat = gsl_matrix_calloc(rows, cols);
sleep(20);
gsl_matrix_free(testmat);
printf("---> Watch eigen matrix memory usage!
");  

Eigen::SparseMatrix<double> mat(rows,cols);         // default is column major
FillMatrix(& mat);
printf("------------------------>DONE
");
return(0);
}

/*-- --*/
int FillMatrix(Eigen::SparseMatrix<double> &mat)
{
int i, j;
Eigen::VectorXd Vres;
mat.reserve(Eigen::VectorXi::Constant(cols,3000));
    for(i=0;i<rows;i++)
    {
        if(i%500==0){printf("i= %i
", i);}
    for(j=0;j<cols;j++)
        {
        if (vee(i,j) != 0){mat.insert(i,j) = vee(i,j);    /*alternative: mat.coeffRef(i,j) += v_ij;*/ }
        }

    }
printf("--->starting compression
");   
mat.makeCompressed();

return(0);
}
/*-- --*/
double vee(int i, int j)
{
    double result = 0.0;

    if(i%4 == 0){result =1.0;}

    return result;
}
/*-- --*/

The error messages are:

ludi@ludi-M17xR4:~/Desktop/tests$ g++ -o eigenfill.x eigenfill.cc -L/usr/local/lib -lgsl -lgslcblas && ./eigenfill.x
eigenfill.cc: In function ‘int main()’:
eigenfill.cc:24:17: error: invalid initialization of non-const reference of type ‘Eigen::SparseMatrix<double>&’ from an rvalue of type ‘Eigen::SparseMatrix<double>*’
 FillMatrix(& mat);
                 ^
eigenfill.cc:12:5: error: in passing argument 1 of ‘int FillMatrix(Eigen::SparseMatrix<double>&)’
 int FillMatrix(Eigen::SparseMatrix<double> & mat);
     ^
ludi@ludi-M17xR

EDIT And it compiles if I write:

FillMatrix(mat);

instead of

FillMatrix(&mat);

I don't understand. Should not the last one be correct?

Answer 1

This is because you are entirely filling one column over four. It looks like the following:

1 0 0 0 1 0 0 0 1 ...
1 0 0 0 1 0 0 0 1 ...
1 0 0 0 1 0 0 0 1 ...
1 0 0 0 1 0 0 0 1 ...
...

Replacing j%4 == 0 by i%4 == 0 will do the job. The memory usage will be exactly: 2500*cols*(sizeof(double)+sizeof(int)) + cols*sizeof(int).

You must also fix the prototype of FillMatrix so that you fill the real matrix and not a temporary:

int FillMatrix(Eigen::SparseMatrix<double> &mat);

Finally, if your goal is really to work with sparse matrices having about 1/4 of non zeros, then a dense representation will likely be faster.