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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Mon, Mar 18, 2013 at 8:15 PM,  <Daniel.Vollmer@xxxxxx> wrote:
> There's one oddity I don't understand: I've had to define
> Eigen::internal::scalar_product_traits for <BlockMatrix, BlockMatrix::Scalar>
> (which I found a bit strange as Matrix<double, 5, 5> * double should already
> be available).

yes, that's to tell that 1) this product is ok, and 2) what is the
return type of this product. This specialization could be provided by
Eigen though.

> A problem I'm currently stuck at is why Eigen doesn't seem to pick up my
> definitions of real() and abs2():
> The calls happen in ConjugateGradient.h:49
>   RealScalar absNew = internal::real(residual.dot(p));
> => residual and p are vectors of my Scalar-type, so vectors of 5x5 matrices for
> example.
> The other culprit is in the next line, i.e. ConjugateGradient.h:50
>   RealScalar rhsNorm2 = rhs.squaredNorm();
> which is from Dot.h:115
>   return internal::real((*this).cwiseAbs2().sum());
> but it seems the cwiseAbs2() again only finds the internal implementation..
>
> The actual compilation errors are as follows (from the internal::real & abs2):
> eigen-3.2.0/Eigen/src/Core/MathFunctions.h:68:12: No viable conversion from 'const BlockMatrixWrapper<Eigen::Matrix<double, 5, 5, 0, 5, 5> >' to 'RealScalar' (aka 'double')
>
> eigen-3.2.0/Eigen/src/Core/MathFunctions.h:263:12: No viable conversion from 'const typename ProductReturnType<Matrix<double, 5, 5, 0, 5, 5>, Matrix<double, 5, 5, 0, 5, 5> >::Type' (aka 'const CoeffBasedProduct<LhsNested, RhsNested, EvalBeforeAssigningBit | EvalBeforeNestingBit>') to 'RealScalar' (aka 'double')
>
> in spite of me having defined
> template<typename B>
> inline typename Eigen::NumTraits<typename B::Scalar>::Real abs2(const BlockMatrixWrapper<B>& x)
> and
> template<typename B>
> inline typename Eigen::NumTraits<typename B::Scalar>::Real real(const BlockMatrixWrapper<B>& x)
>
> Attached is my current code (which doesn't compile due to the mentioned problems).

I'm afraid this won't suffice to make it work. For instance if you
look at the CG algorithm, you see it needs some intermediate vectors.
If we declare them as the right-hand-side, in your case a vector of
block, then expressions like "Scalar alpha = absNew / p.dot(tmp);"
cannot work because   p.dot(tmp) must return a double not a block.
Here "Scalar" should also be a double and not a block. Now if we
declare the temporary vectors as vectors of doubles, then we have
types and dimensions mismatch.

What you really need is a SparseBlockMatrix wrapper whose Scalar type
is double, rows(), cols() returns the dimension at the scalar level
(not the number of blocks, etc.). Then you'll have to overload the
features we want to support for them to cast them in terms for
(sparse) matrix of blocks. For the CG, overloading operator* with a
dense RHS should be enough.

gael


>
>
> Best regards
>
> 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