Re: [eigen] Way to obtain AdjointReturnType?

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi Rhys,

I didn't look into how to fix your errors, but I'd like to say the
following so you check if you really need to bother with all this,
first.

In the development version, the "views" like SelfAdjointView are not
sub-classes of MatrixBase.

This means that SelfAdjointView doesn't have the generic matrix API.
It only has a little selfadjoint-specific API.

So, I didn't really expect that people would write functions returning
a SelfAdjointView. Even if the returned matrix happens to be
selfadjoint. Because that would seriously limit what can be done with
that object. It would make sense if you wanted to only fill one
triangular part and let the other triangle be implicit from
selfadjointness, but apparently that's not what you're doing.

Instead, I'd expect such a function to return a plain matrix, and then
when you want to use it in an operation where its selfadjointness
matters, do

matrix.selfAdjointView().someOperation()

In other words, SelfAdjointView is meant as a type of temporary
wrapper objects used to call the selfadjoint-specific API.

Does this work for you?

Benoit

2009/9/7 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
>
>