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?

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 .

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.

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";
}


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

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.

Karma can be awarded by following the Karma link in the postbit or on the user's profile, then selecting "Rate this user".

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.

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?

Or, in a similar question, can I write a function which will happily accept a true Eigen::MatrixXf* OR the CWisenullaryop<Eigen::Matrixxf>* object?

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!

For the HilbertOp case above, what would be an example definition of the SIMD version of the function? Thanks in advance.

alvaros wrote:For the HilbertOp case above, what would be an example definition of the SIMD version of the function? Thanks in advance.

bump