Re: [eigen] Initializing a sparse matrix - surprisingly slow

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

Many thanks Gael!

setFromTriplets is very comfortable and only a bit slower then using
insert + A.reserve(Eigen::VectorXi::Constant(n,5));

But, it uses 4 times as much memory:
VectorXi::Constant(n,5)  uses  n integers  =  4n  bytes
n Triplets use  n x (2*sizeof(int)+sizeof(double)) = 16n bytes.

(in my easy case the elements of the matrix can only have  one of 4  
values.)

Perhaps there should be a warning when using "insert" without the  
"correct"
"reserve". On a matrix of dimension a million it's slower by a factor
of more than 10000 !

Helmut.

On 10/02/2012 01:28:22 PM, Gael Guennebaud wrote:
> Hi,
>
> in order to not bother about matrix assembly, the recommended way is
> to use a triplet list as explained in the tutorial:
>
> http://eigen.tuxfamily.org/dox-devel/TutorialSparse.html#TutorialSparseFilling
>
> Currently, .insert() always transform the matrix into a non-compressed
> form by calling
> this->reserve(Eigen::VectorXi::Constant(this->outerSize(),2)); that
> makes insertions extremely slows fo ran obvious reasons. In your case,
> to fully take advantage of a fully coherent insertion strategy you
> would have to use the low level API which is not documented anymore
> because I prefer to see standard users use the setFromTripplets API:
>
>  for (int j= 0; j < n; j++) {
>     A.startVec(j);
>     A.insertBack(i_1,j) = ...;
>     A.insertBack(i_2,j) = ...;
>     ...
> }
> A.finalize();
>
> with the condition that: i_k < i_{k+1}.
>
> gael
>
> On Tue, Oct 2, 2012 at 11:53 AM, Helmut Jarausch
> <jarausch@xxxxxxxxxxxxxxxxxxx> wrote:
> > Hi,
> >
> > probably I haven't understood the docs:
> >
> > #include <Eigen/Sparse>
> > using Eigen::SparseMatrix;
> > typedef SparseMatrix<double>  SpMatrix;
> >
> > #define Slower_by_a_factor_of_1000
> >
> > SpMatrix Laplace( int m ) {  // 2D Laplace : 5 elements per  
> row/column at
> > most
> >   const int  n = m*m,  n0= n-1, n_m = n-m,
> >              NumNonZeros= 5*n;  // number of non-zeros
> >   SpMatrix A(n,n);
> >
> > #ifdef Slower_by_a_factor_of_1000
> >   A.reserve(NumNonZeros);  // 1000 times slower for n=9801 <<<<<  
> WHY ???
> > #else
> >   A.reserve(Eigen::VectorXi::Constant(n,5));
> > #endif
> >
> >   double d = 4,  b = -1, c = -1;
> >
> >   for (int j= 0; j < n; j++) {
> >     if ( j >= m ) A.insert(j-m,j)= c;
> >     if ( j > 0 &&  j%m != 0 ) A.insert(j-1,j)= b;
> >     A.insert(j,j)= d;
> >     if ( j < n0 && (j+1)%m != 0 ) A.insert(j+1,j)= b;
> >     if ( j < n_m ) A.insert(j+m,j)= c;
> >   }
> >   return A;
> > }
> >
> > This is with the GIT-version of Eigen.
> >
> > Many thanks for an explanation,
> > Helmut.
> >
> >