Re: [eigen] Way to obtain AdjointReturnType?

[Rhys Ulerich <rhys.ulerich@xxxxxxxxx>]

Hi Gaël,

> I guess you are using the 2.0 branch then you probably want to return a:
>
> Flagged<AnEigenMatrixType,SelfAdjointBit,0>
> ...
> Part<AnEigenMatrixType,SelfAdjoint | UpperTriangular> // or: | LowerTriangular

> ...
> SelfAdjointView<AnEigenMatrixType,UpperTriangular> // or: | LowerTriangular
> ...
> if you want more information please be more specific about your use case.

I'm trying to use the SelfAdjointView and Flagged approaches that you suggested without much luck.  I'd appreciate any hints you could give about how to express what I'm after.  The more-specific use case is computing the viscous stress tensor in a compressible Navier-Stokes simulation:

template<typename Scalar>
Eigen::SelfAdjointView<Eigen::Matrix<Scalar,3,3>,Eigen::UpperTriangular> tau(
        const Scalar &mu,
        const Scalar &lambda,
        const Scalar &div_u,
        const Eigen::Matrix<Scalar,3,3> &grad_u)
{
    return mu*(grad_u + grad_u.transpose())
        + lambda*div_u*Eigen::Matrix<Scalar,3,3>::Identity();
}

My template parameter will always be a real floating point type, and so the returned matrix will always be symmetric.  The guaranteed symmetry of the returned, real-valued matrix is what I'm trying to capture using the SelfAdjointView.

When I instantiate the above template using Eigen 2.0.5 with the Intel 10.1 compiler I see

../../sz/suzerain/classical.hpp(327): error: no suitable user-defined conversion from
         "const Eigen::CwiseBinaryOp<Eigen::ei_scalar_sum_op<double>, Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<double>, Eigen::CwiseBinaryOp<Eigen::ei_scalar_sum_op<double>, Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Transpose<Eigen::Matrix<double, 3, 3, 0, 3, 3>>>>, Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<double>, Eigen::CwiseNullaryOp<Eigen::ei_scalar_identity_op<double>, Eigen::Matrix<double, 3, 3, 0, 3, 3>>>>" to
         "Eigen::SelfAdjointView<Eigen::Matrix<double, 3, 3, 0, 3, 3>, 1024U>" exists
     return mu*(grad_u + grad_u.transpose())
            ^
         detected during instantiation of "Eigen::SelfAdjointView<Eigen::Matrix<Scalar, 3, 3, 0, 3, 3>, 1024U> pecos::suzerain::orthonormal::rhome::tau(const Scalar &, const Scalar &, const Scalar &, const Eigen::Matrix<Scalar, 3, 3, 0, 3, 3> &) [with Scalar=double]" at line 345 of "../../sz/tests/test_classical.cpp"

My second attempt at using SelfAdjointView was...

template<typename Scalar>
Eigen::SelfAdjointView<Eigen::Matrix<Scalar,3,3>,Eigen::LowerTriangular> tau(
        const Scalar &mu,
        const Scalar &lambda,
        const Scalar &div_u,
        const Eigen::Matrix<Scalar,3,3> &grad_u)
{
    using namespace Eigen;

    const SelfAdjointView<Matrix<Scalar,3,3>,LowerTriangular> tmp(
                grad_u + grad_u.transpose());

    return mu*tmp + lambda*div_u*Eigen::Matrix<Scalar,3,3>::Identity();
}

... and it wasn't any more successful:

.../../sz/suzerain/classical.hpp(332): error: no operator "*" matches these operands
            operand types are: const double * const Eigen::SelfAdjointView<Eigen::Matrix<double, 3, 3, 0, 3, 3>, 2048U>
      return mu*tmp + lambda*div_u*Eigen::Matrix<Scalar,3,3>::Identity();
               ^
          detected during instantiation of "Eigen::SelfAdjointView<Eigen::Matrix<Scalar, 3, 3, 0, 3, 3>, 2048U> pecos::suzerain::orthonormal::rhome::tau(const Scalar &, const Scalar &, const Scalar &, const Eigen::Matrix<Scalar, 3, 3, 0, 3, 3> &) [with Scalar=double]" at line 345 of "../../sz/tests/test_classical.cpp"

I've not had any better luck on 2.0.5 using Flagged.  The templated return type

Eigen::Flagged<Eigen::Matrix<Scalar,3,3>,Eigen::SelfAdjointBit,0> tau(/*...*/)

results in

../../sz/suzerain/classical.hpp(327): error: no suitable user-defined conversion from
          "const Eigen::CwiseBinaryOp<Eigen::ei_scalar_sum_op<double>, Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<double>, Eigen::CwiseBinaryOp<Eigen::ei_scalar_sum_op<double>, Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Transpose<Eigen::Matrix<double, 3, 3, 0, 3, 3>>>>, Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<double>, Eigen::CwiseNullaryOp<Eigen::ei_scalar_identity_op<double>, Eigen::Matrix<double, 3, 3, 0, 3, 3>>>>" to
          "Eigen::Flagged<Eigen::Matrix<double, 3, 3, 0, 3, 3>, 512U, 0U>" exists
      return mu*(grad_u + grad_u.transpose())
             ^

I haven't been able to get your Part suggestion going either.

Thank you for your time,
Rhys

> > I've got a templated function with a signature like
> >    ReturnType foo(AnEigenMatrixType bar)
> > where AnEigenMatrixType can be an instantiation of Eigen::Matrix<...>.  I'd

> > like ReturnType to be AnEigenMatrixType with its SelfAdjointBit set because
> > I know the computation in foo() guarentees me a self-adjoint result.
>
> I guess you are using the 2.0 branch then you probably want to return a:

>
> Flagged<AnEigenMatrixType,SelfAdjointBit,0>
>
> or
>
> Part<AnEigenMatrixType,SelfAdjoint | UpperTriangular> // or: | LowerTriangular
>
> according to want you want to achieve. Basically the former as no

> effect in Eigen while the other tells only the upper or lower
> triangular part is meaningful.
>
> However the support for selfadjoint matrices in the 2.0 branch is not
> very good, and in the devel branch the self adjoint bit flag, .part()

> and Flagged are now deprecated in favor of a SelfAdjointView wrapper
> class. Then your return type would be:
>
> SelfAdjointView<AnEigenMatrixType,UpperTriangular> // or: | LowerTriangular

>
> if you want more information please be more specific about your use case.
>
> cheers,
> gael.
>
> --
> Gaël Guennebaud
> Iparla - INRIA Bordeaux
> (+33)5 40 00 37 95

> > I noticed that MatrixBase<Derived> has AdjointReturnType used in adjoint(),
> > but my attempts like Eigen::MatrixBase<Eigen::Matrix<Scalar,3,1>
> >>::AdjointReturnType haven't worked (public-vs-private, but it could be me

> > being dumb with templates).  Functionality like part<Eigen::SelfAdjoint>
> > doesn't quite do what I'd like because I want to propagate the
> > self-adjointedness of ReturnType back out to other routines.

> >
> > What am I missing?
> >
> > Thanks,
> > Rhys