Hello,
I am replacing a no-longer-supported Linear Algebra library (RogueWave Math.h++) with Eigen in a large simulation. There are some matrix and vector functions that I require that aren't in Eigen, such as abs() and sort(). These functions are called with varying argument types (e.g. VectorXd, VectorXi, MatrixXd, etc.), so I need a templatized function to handle all cases.

An example with the method I attempted is:

Code: Select all

#include <Eigen\Core>
   USING_PART_OF_NAMESPACE_EIGEN

// Example templatized function
template<typename Derived>
Derived abs( const Derived& v )
{   
   Derived vabs = v;
   for(int i=0; i<vabs.size(); i++)
      if(vabs(i)<0)
         vabs(i) = -vabs(i);
   return vabs;
};

int main(void)
{
   VectorXd x1(3); x1.fill(-1.);
   VectorXd x2(3); x2.fill(-2.);
   MatrixXd A(3,3); A.fill(-3.);

   VectorXd y;
   MatrixXd B;

   // Templatized function works when type is expressly stated (i.e. <type>)
   y = abs<VectorXd>(x1); // Works.
   B = abs<MatrixXd>(A); // Works.
   y = abs<VectorXd>(x1-x2); // Works.
   y = abs<VectorXd>(A*x1); // Works.

   // Templatized function can fail...
   y = abs(x1); // Works. (assumes VectorXf at compile)
   B = abs(A); // Works. (assumes MatrixXf at compile)
   y = abs(x1-x2); // Fails. (assumes CwiseBinaryOp<...> at compile)
   y = abs(A*x1); // Fails. (assumes Product<...> at compile)

   return 0;
}

The above templatized abs() works if I explicitly state the type. But it does not handle "x1-x2" and "A*x" as arguments when the type isn't explicitly stated.

I'd be happy to use the above and explicitly state the type (e.g. abs<VectorXf>(x1-x2) ), if abs(x1-x2) yielded a compile-time or run-time error. The problem is, it doesn't. Instead, the program just stops. (I believe it gets a stack overflow due to infinite recursion, as the CwiseBinarOp<> doesn't support calling vabs(i) in the function abs().) Eventually, one of our developers will write "abs(x1-x2)", and the simulation will just stop in a manner that will be very difficult to debug.

Any suggestions this community has will be very appreciated!

As you realized, the expression A*x1 has type Product<...>. Your function should return an object of type VectorXd. This type is available as the member PlainMatrixType. Here is a version of your abs function with minimal changes that should work.

Code: Select all

template<typename Derived>
typename Derived::PlainMatrixType abs( const Derived& v )
{   
   typename Derived::PlainMatrixType vabs = v;
   for(int i=0; i<vabs.size(); i++)
      if(vabs(i)<0)
         vabs(i) = -vabs(i);
   return vabs;
};


Incidentally, the absolute value exists in Eigen 2.0 as Cwise::abs(); see http://eigen.tuxfamily.org/dox/classEigen_1_1Cwise.html#abf73276fc968045bd7d6b9dbb1fe3b72
So you can write

Code: Select all

y = (x1-x2).cwise().abs();


You may need to "#include <Eigen\Array>".

Sweet! PlainMatrixType worked perfectly and is exactly what I needed to know. And also good to know about the additional cwise() functionality that I had overlooked. Many, many thanks!