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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

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.
>
>