Re: [eigen] SparseMatrix with non POD-"Scalar" type (e.g. SparseBlockMatrix) yields unnecessary initialisations of Scalar-type?

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

What you describe is exactly what I have in mind because it's what I
need too, but of course I'd also like to have the possibility to have
all blocks the same fixed size.

gael

On Sun, Mar 17, 2013 at 9:25 PM, Sameer Agarwal
<sameeragarwal@xxxxxxxxxx> wrote:
> Gael,
> Sparse block matrices in Eigen would be a really great addition.
>
> Is the plan to have all blocks of the same size?
>
> One direction of interest to us would be sparse matrices that support
> row and column block size vectors, e.g
>
> consider "3x3" sparse block matrix
>
> col_blocks = [2 4 2]
> row_blocks = [3 4 5]
>
> Then the block structure will look like
>
> M = [3x2 3x4 3x2]
>        [4x2 4x4 4x2]
>        [5x2 5x4 5x2]
>
> Sameer
>
>
> On Sun, Mar 17, 2013 at 11:38 AM, Gael Guennebaud
> <gael.guennebaud@xxxxxxxxx> wrote:
>> Hi Daniel,
>>
>> On Sat, Mar 16, 2013 at 7:00 PM,  <Daniel.Vollmer@xxxxxx> wrote:
>>
>>> I've attached my experiment that essentially uses a dense 5x5 matrix as the "Scalar" type. Basically, this works (which is already pretty nice) but I have noticed two things which seem sub-optimal, and so I was curious whether I'm missing something or doing it wrong.
>>
>> adding an explicit support for sparse block matrices is something we
>> planed to address soon. Your example is nice starting
>> proof-of-concept. To get it right, we certainly need a new
>> SparseBlockMatrix class and some way to specify which algorithm to use
>> for a/b, i.e., a dense Cholesky in SimplicialCholesky, LU in SparseLU,
>> etc.
>>
>>> 1) The Scalar-type used in SparseMatrix seems to need an abs() function in order to (I assume) only iterate over non-zero elements (AmbiVector.h:308 & 335). Oddly enough, my abs-version was needed to compile, but never got called.
>>> It might be advantageous to provide an iterator that iterates over all elements without making such a guarantee, but then the point shouldn't matter once it is compressed.
>>
>> abs is probably not the only free fonction you/we will have to
>> specify, but also operator/, and in some cases sqrt (LLT). In this
>> precise case of AmbiVector, is only needed for non-conservative sparse
>> matrix products, i.e.:  (A*B).prune(eps);
>>
>>> 2) Many times the single-argument constructor is called for the Scalar-type of SparseMatrix in instances where I don't think this is necessary. As this becomes more and more expensive with a more heavy-weight Scalar type, I'd like to avoid this where possible, for example
>>> - SparseMatrix.h:1108 (insertUncompressed): Element is initialised to 0 even though it is later overwritten by the inserted value. Worked fine for me if I removed the = 0.
>>> - SparseMatrix.h:382 (insertBackByOuterInner) and SparseMatrix.h:392 (insertBackByOuterInnerUnordered): Element is initialised to 0 even though it is later overwritten by the inserted value. I hacked this by making a variant of CompressedStorage::append() that only takes an Index and not a value-ref and calling that instead (because again, the value is later overwritten via the reference we return)
>>> - AmbiVector.h:197/206/239 (coeffRef): I don't quite understand this, but I'm not sure the element needs to be initialised here either.
>>
>> We need these initializations because it is allowed to do:
>>
>> mat.insert(i,j);
>>
>> to insert an explicit zero. One solution could be to return a wrapper
>> object that would set the wrapped value to zero in the destructor if
>> operator= has not been called, but that sounds a bit overkill.
>>
>>> - AmbiVector.h:171 (setZero): Shouldn't / couldn't this use assignment from m_zero instead?
>>
>>
>> Actually, all these constant-to-Scalar conversions are problematic for
>> any heavy-weight scalar types (AutoDiff, multi-precision, etc.), and
>> an old idea was to add in NumTraits a Constant (maybe that name is not
>> good) type for that usage. This way if custom scalar type is fully
>> compatible with double, then NumTraits<Custom>::Constant could be
>> double. For fixed-size matrices it could be a CwiseNullaryOp<> thus
>> making the assignment more efficient and optimizable by the compiler.
>>
>>> Is this approach (using a custom scalar type) still the recommended way to create sparse block-matrices?
>>
>> As I said, there is currently no recommended way to create
>> sparse-block matrices but this is something we'd like to have, and
>> we'll probably work on it soon. So if you're willing to share your
>> advances in this direction and participate in any way to this
>> development, you're very welcome!!
>>
>>> Can we reduce the overhead of these initialisations or make them optional?
>>
>> Do you have some measure of the overhead?
>>
>> gael
>>
>>>
>>> I'd welcome any feedback. :)
>>>
>>> Thanks!
>>>     Daniel Vollmer
>>>
>>> --------------------------
>>> Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR)
>>> German Aerospace Center
>>> Institute of Aerodynamics and Flow Technology | Lilienthalplatz 7 | 38108 Braunschweig | Germany
>>
>>
>
>