Re: [eigen] Sparse usage feedback

[<bernard.hugueney@xxxxxxxxxx>]

Hi Geal,



Thanks for your answers. I reply "inline".



> 

> I could consider to give it a try if and only if the begin()/end()

> paradigm would bring any additional benefits than just being closer to

> the STL. The current approach has the advantage to be simpler to use

> (e.g., no need to declare an end iterator). Well, actually I'm

> thinking that having an access to the end() might be useful but I'm

> still a bit worried about the performance and the API verbosity.

> 



In fact, my wish for more idiomatic iterators doe not stems from the

abstract love of iterator pairs, but from the practical desire to reuse

idiomatic code (such as the STL itself) on sparse matrix data ! :)

[more below ]



> Here is an example of what could be the API:

> 

> SparseMatrix m;

> 

> for(int j=0;j<m.outerSize();++j)

> {

>   SparseMatrixType::InnerIterator end(m.innerVector(j).end());

>   for(SparseMatrixType::InnerIterator it(m.innerVector(j).begin());

> it!=end; ++it)

>     ...;

> }



That would be good. I also pondered about an API inspired from

tr1::unordered_map :

http://publib.boulder.ibm.com/infocenter/comphelp/v9v111/index.jsp?topic=/com.ibm.xlcpp9.aix.doc/standlib/stl_unordered_map.htm



with .begin(IndexType) & .end(IndexType).



[more below]

> 

>> 2°) The iterator do not expose row() and col() through an iterator

>>

>> interface

>>

>> (i.e. with operator*()) which make it impossible to manipulate

sequences

>>

>> of

>>

>> coords with the standard algorithms operating on iterator sequences.

> 

> Sorry, but I don't really understand what you are looking for here,

> can you be more specific please?



I meant that you can have access to all the relevent info with your

current API,

but not through an "itarator API" because you get them by calling member

functions 

it.value();

it.row();   // row index

it.col();   // col index (here it is equal to k)

it.index(); // inner index, here it is equal to it.row()



which means that you cannot replace your loops with calls to std

algorithms.

I.e, I'd like to rewrite your previous example as :

SparseMatrix m;

 

for(int j=0;j<m.outerSize();++j)

  std::for_each(m.innerVector(j).begin(), m.innerVector(j).end(), …);



But it means that the iterators must return the relevent info through the 

idiomatic operator*() instead of custom member functions (.value(),…).



I wrote iteratorS because as you only have one operator*(), you'd have to

give more than one iterator, or you'd make your user pay for unused data.

(i.e you could give rowsCols(j).begin()/.end() with operator*() returning

and std::pair<IndexType, IndexType>

and values(j).begin()/.end() returning a Scalar.



Failing to do that means giving up on the std algorithms (or having to

rewrite them

as does GIL : 

http://stlab.adobe.com/gil/html/giltutorial.html#STLEquivalentsSec



Btw, I'll also investigate whether it would be prohibitively expensive to

allow

 std::for_each(m.begin(), m.end(), …);

to walk the whole matrix.(IMVHO, for sparse matrices, you don't usually

fit many 

matrices in L1 cache so the extract checks will be lost in the data fecth

latency).





> 

>> 3°) This is not possible to construct a Sparse from a Dense

>>

> 

> This is already possible:

> 

> sparse = dense.sparseView();

> sparse = dense.sparseView(ref,epsilon); // skip coefficients smaller

> than ref*epsilon



Thanks. I missed (and I think it should deserve some ref in 

http://eigen.tuxfamily.org/dox-devel/TutorialSparse.html



However, I think that Sparse classes should have constructors taking dense

matrices

calling sparseview() because if you can construct an object of class A

from 

an object of classe B, the idiomatic way to do it in C++ is to have a A(B

const&)

constructor.





> 

> 

>> 4°) I think some scalar operation traits are missing

>>

>> (i.e. trying to multiply an int matrix with a double fails to find

>>

>> ‘struct Eigen::ei_scalar_product_traits<double, int>’ )

> 

> This is because this is explicitly not allowed. You have to explicitly

> cast the int matrix to double. Our opinion is that such an automagic

> cast would be dangerous.



When sizeof(int)*8 <= std::numeric_limits<double>::digits I believe that

such

a cast is not dangerous (and I have 32 bits for ints and 53 bits for

double).



However, I think maybe a catch all template could be used to give an hint

to the

problem with an helpful message in the source code of this template

(designed to trigger an error)

such as /* no implicit cast provided, use the .cast<> template member

function 

to mix those types*/).



> 

>> 5°) I think that the typedefs in i.e.  SparseLU<> should be public

>>

>> I may need to know the LMatrixType & cie to make some metaprogramming

of

>>

>> my own

>>

>> (or did I missed some traits struct ?)

> 

> Yes you are right. All these decompositions/solvers require a lot of

> polishing, and actually I've already started a rewrite of them with a

> slightly simpler/safer API. The main change is that I'm removing these

> Backend template parameters.



I am very interested in this !

Where is it availbale ?



Currently, I'm giving the backend choice at runtime and deriving

the implementations with dynamic dispatching.

    template<typename Backend> struct LUFact : LUFactBase



and I convert the returned matrices types to a Backend-independant

sparse type.



As I wrote above, I think that sparses matrices (contrary to dense

matrices)

are never small enough to mandate the all same optimizations as dense

matrices.

For instance, I can't think of a relevent use case forbidding dynamic

dispatching

(maybe using covariant return type if various backend return different

concrete

implementations for the result matrices of the LU decompositions).



Best regards,



Bernard