Reply to topic

wrap a function as an eigen matrix

johnm1019
Registered Member
Posts
45
Karma
0

wrap a function as an eigen matrix

Thu Jul 26, 2012 9:01 am
Hi, we currently generate our system matrix in memory (as Eigen::Matrix), then run our solver.

In the game of speed/memory tradeoff's we want to generate a version of our code that computes the system matrix on the fly. We can wrap this as a function which looks like sysMat(i,j). Ideally though, we would also wrap it as an Eigen object so that it fits easily into our existing solver and more importantly takes advantage of SSE optimizations.

Is it possible to wrap an Eigen::Matrix around a function which takes the indecies you want as arguments?
jitseniesen
Registered Member
Posts
185
Karma
2
This is perhaps best done using CwiseNullaryOp. Unfortunately, the documentation for this part of Eigen is quite lacking. You can see some examples in ticket 457 at http://eigen.tuxfamily.org/bz/show_bug.cgi?id=457 .
johnm1019
Registered Member
Posts
45
Karma
0
Thanks for that. From the comments it appears to do what I want, but the specific declaration of the header I need I cannot determine.

I want a function which looks like an Eigen::MatrixXf, and only for float's, and always 2D layout.
jitseniesen
Registered Member
Posts
185
Karma
2
Here is an example which you can hopefully adapt to your use case. It computes the Hilbert matrix (see http://en.wikipedia.org/wiki/Hilbert_matrix):
Code: Select all
#include <iostream>
#include "Eigen/Core"

struct HilbertOp
{
  // This is the function that computes the matrix entries on the fly
  template<typename Index>
  const float operator() (Index row, Index col) const
  {
    return 1 / (float(row+1) + float(col+1));
  }
};

namespace Eigen {
  namespace internal {
   
    // This lists various properties of the HilbertOp functor:
    // Cost = estimate of cost of calling the functor
    // IsRepeatable = whether functor gives same result if called
    //                multiple times with the same parameters
    // PacketAccess = whether vectorized variant of functor is implemented
    template <> struct functor_traits<HilbertOp> {
      enum { Cost = 10,
             IsRepeatable = true,
             PacketAccess = false
      };
    };

    // This specifies that the functor is to be called with two parameters
    // (row index and column index)
    template<> struct functor_has_linear_access<HilbertOp> {
      enum { ret = 0 };
    };
  }
}

int main()
{
  // This creates a 5-by-5 matrix. The entries are not yet computed.
  Eigen::CwiseNullaryOp<HilbertOp, Eigen::MatrixXf> mat(5, 5, HilbertOp());

  // Print the matrix. This computes the entries by calling HilbertOp::operator()
  std::cout << mat << "\n";
}
johnm1019
Registered Member
Posts
45
Karma
0
This is awesome, thank you! Eigen is amazing.

Can you elaborate on how to set the Cost parameter, as well as in what situations supporting PacketAccess is advisable?
There is no way to hand-optimize my function such that calling myFunc(1,i) for i=0:3 would be any faster than calling myFunc(1,0,1,2,3) -- so I assume this would be the case where PacketAccess is used? (if I could hand code a 4-wide version that was faster than calling 1-wide 4x?).

Thanks

PS - I want to give you Karma but I don't see how I can do that :-\
User avatar ggael
Moderator
Posts
2193
Karma
15
OS
Cost permits to control unrolling and intermediate evaluations. The cost of one product or or one addition of floats or doubles is 1. If your function is quite complex, just pust an arbitrary high number like 1000.

PacketAccess means you can provide an SIMD version.
User avatar bcooksley
Administrator
Posts
18556
Karma
83
OS
Karma can be awarded by following the Karma link in the postbit or on the user's profile, then selecting "Rate this user".


System Settings and Device Actions KCM maintainer
Image
johnm1019
Registered Member
Posts
45
Karma
0
bcooksley wrote:Karma can be awarded by following the Karma link in the postbit or on the user's profile, then selecting "Rate this user".

Thanks for the tip!
I must say - this is outrageously difficult to find using the default theme.
alvaros
Registered Member
Posts
2
Karma
0
Hi, I have done the above mentioned and have already wrapped my coefficient generator function with an Eigen::MatrixXf using the CwiseNullaryOp, works awesome. Now I can compute any coefficient at will by calling sysMat(i,j).

What I need to do now, is to return my CwiseNullaryOp<coeffxy, Eigen::MatrixXf>* mat as an Eigen::MatrixXf* and use this return value as an argument for the solver:

Code: Select all
// system matrix containing the cwisenullaryop
CwiseNullaryOp<coeffxy, Eigen::MatrixXf>*  SystemMatrix::myNullaryOp;
Eigen::MatrixXf* SystemMatrix::getMatrix()
{
    return  myNullaryOp;
};

// solver definition
Solver::solve(Eigen::MatrixXf* matrix){..};

// solver call
theSolver.solve( sysMat.getMatrix() );


In essence, is it possible to make a Eigen::CwiseNullaryOp look exactly like an Eigen::MatrixXf ? This of course without allocating the whole matrix but rather keeping the otf capabilities of the nullaryOp?
johnm1019
Registered Member
Posts
45
Karma
0
Or, in a similar question, can I write a function which will happily accept a true Eigen::MatrixXf* OR the CWisenullaryop<Eigen::Matrixxf>* object?
User avatar ggael
Moderator
Posts
2193
Karma
15
OS
You need to template your function. See:

http://eigen.tuxfamily.org/dox/TopicFun ... Types.html

Otherwise, the coeff(i,j) method would have to be virtual, that would lead to extremely slow computations!
alvaros
Registered Member
Posts
2
Karma
0
For the HilbertOp case above, what would be an example definition of the SIMD version of the function? Thanks in advance.
johnm1019
Registered Member
Posts
45
Karma
0
alvaros wrote:For the HilbertOp case above, what would be an example definition of the SIMD version of the function? Thanks in advance.

bump

 
Reply to topic

Bookmarks



Who is online

Registered users: AGB, Azzart, Baidu [Spider], Bing [Bot], Bingy, colomar, Exabot [Bot], garthecho, Google [Bot], google01103, GreatEmerald, jensreuterberg, jstaniek, kainz.a, koriun, Kver, La Ninje, lazyit, Majestic-12 [Bot], MSN [Bot], Nuc!eoN, psbot [Picsearch], rulet111, scummos, SecretCode, Tuukka, undoshmulhe, Uri_Herrera, Yahoo [Bot]