Re: [eigen] Clamping the coeffs of a matrix

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

2010/1/23 Alexandre Devert <marmakoide@xxxxxxxx>:
> Ho hai,
>
>  I'm working on reaction-diffusion problems. I need to clamp the coefficients of a dense matrix, as following :
> W(i, j) = max(a, W(i, j))
>
>  A more generalization would be
> W(i, j) = min(b, max(a, W(i, j)))  for clamping in [a ,b]
>
> Questions :
>  1) Can I do it in a better way than a naive loop over the coefficients ?

Yes, using a custom functor with unaryExpr(). By chance, this is
exactly what our example for that does:
http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#a1d1835f182c0141dc073a1045c8a2e9e

This will allow to traverse the array only once, and unroll if appropriate.

>  2) I know that the max and min operator are in the SSE instruction set, thus, can I benefit from it ?

Yes.
First of all, notice that in Eigen/src/Core/arch/SSE/PacketMath.h we
have wrapper functions ei_pmin() and ei_pmax() that wrap the SSE
intrinsic functions for that in a portable way.
Now you want to make use of them in your custom functor. To do that,
just copy what we do in other functors that are already vectorized.
See for example in Eigen/src/Core/Functors.h the ei_scalar_sum_op:

template<typename Scalar> struct ei_scalar_sum_op EIGEN_EMPTY_STRUCT {
  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const
Scalar& b) const { return a + b; }
  template<typename PacketScalar>
  EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar&
a, const PacketScalar& b) const
  { return ei_padd(a,b); }
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_sum_op<Scalar> > {
  enum {
    Cost = NumTraits<Scalar>::AddCost,
    PacketAccess = ei_packet_traits<Scalar>::size>1
  };
};

As you can see:
 - you must implement the packetOp() method
 - you must announce that your functor supports SIMD by setting
PacketAccess to 1 in a specialization of ei_functor_traits.
    (notice that here ei_packet_traits<Scalar> is the number of
Scalar's in a 128bit packet)

Cheers,
Benoit