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

[<Daniel.Vollmer@xxxxxx>]

Hello Gael,

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

I have tried a different approach (see attached files once more), but I
am ignoring any solvers for the moment.

Essentially, I declare
	namespace internal
	{
		template<typename BMA, typename BMB>
		struct scalar_product_traits<BlockMatrixWrapper<BMA>, BlockMatrixWrapper<BMB> >
		{
			typedef BlockMatrixWrapper<Eigen::Matrix<typename scalar_product_traits<typename BMA::Scalar, typename BMB::Scalar>::ReturnType, BMA::RowsAtCompileTime, BMB::ColsAtCompileTime> > ReturnType;
		};
	};
This covers both Matrix * Matrix as well as Matrix * Vector products and makes sure
the resulting block size of a product is correct. The only difficulty with this
(and the single change needed in Eigen to get this to work) was declaration of alpha
in SparseDenseProduct.h (and probably in equivalent places for dense * dense
and sparse * sparse):
  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
  {
    for(Index c=0; c<rhs.cols(); ++c)
    {
      for(Index j=0; j<lhs.outerSize(); ++j)
      {
        typename Res::Scalar rhs_j = /*alpha **/ rhs.coeff(j,c);
        for(LhsInnerIterator it(lhs,j); it ;++it)
          res.coeffRef(it.index(),c) += it.value() * rhs_j;
      }
    }
  }
alpha is declared as Res::Scalar (which in my example is Eigen::Matrix<double, 5, 1>)
and rhs.coeff() is also Eigen::Matrix<double, 5, 1> so the product doesn't make sense.
If we could make sure this alpha is a REAL scalar (I just commented it out as you can see),
then this case would work quite nicely.


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

Attachment: BlockMatrixWrapper.h
Description: BlockMatrixWrapper.h

Attachment: main.cpp
Description: main.cpp