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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Tue, Oct 2, 2012 at 1:50 PM, Helmut Jarausch
<jarausch@xxxxxxxxxxxxxxxxxxx> wrote:
> 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.

Not exactly:
* VectorXi::Constant(n,5) cost 1 integer,
reserve(VectorXi::Constant(n,5) ) allocate (n + 5*n) integers + 5*n
double =>64n Bytes
* A tripplets of 5*n entries cost 80n Bytes + the matrix itself 64n B
=> a total of ~144n Bytes. Once the matrix is assembled, the tripplet
list can be freed, and both strategies become even.

>
> (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 !

I guess we could easily come up with a better heuristic in the case
the reserve is inappropriate. The only problem with complex heuristic
is that it's difficult to predict the behavior.

gael

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