specialise for loop with visitor pattern?

Question

Let's say that I have a struct named like this one:

struct OneBodyMatrixElements {
  static const bool is_selfadjoint = 1;
  static const bool property2 = 0; 
  // ........

  static value_type val( const std::size_t i, const std::size_t j ) {
    value_type ret = 0.;
    if( i == j )  ret = 1.;
    if( i == j + 1 || i == j - 1 ) ret = -1.;
    if( i == 0 && j == 6 ) ret = -1.;
    if( i == 6 && j == 0 ) ret = -1.;
    return static_cast<value_type>(ret);
  }
};

The static variable define different properties of this class and the static function val return the value depending on (i,j). This actually represents a matrix of size 7 whose scalar type is value_type (either float, double, std::complex).

In another part of my code, I need to loop over (i,j) for i,j=0,...,N (here for example N=7):

for( std::size_t i = 0; i < N; i++ ) {
  for( std::size_t j = 0; j < N; j++ ) {
    if( "OneBodyMatrixElements::val( i, j ) is non zero" ) {
      //...
      OneBodyMatrixElements::val( i, j ) * func( i, j );
      //...
    }
}

I have an if close to check if the OneBodyMatrixElements is zero because the multiplication just below is expensive and does not need to be computed if the prefactor is zero.

My question is the following: Since I know the properties of OneBodyMatrixElements, i.e. selfadjoint or diagonal, I also know the location of zero elements, how can I optimise the two for loops to iterate only over non zero elements and also take into account the selfadjoint property?

The obvious answer to this question would be to use Eigen SparseMatrix but I is not what I want because

• I do not want to actually store the matrix
• I actually have another struct with four indices (i,j,k,l) and loop over these four indices. I therefore need something generic that also work for four indices. In Eigen, there is no support for SparseTensors.

I did have a look to the visitor pattern but I'm not sure that this is what is need. Any comments on the code is also welcome and a solution which actually store the values in the memory in a sparse form is acceptable. A solution involving only iterations over non-zero elements is also acceptable. Thank you.

EDIT: Answer:

Thank you very much for your answers. I did implement a custom iterator as explained below. For the sake of completeness, here is my version:

class tridiag_iterator {
  std::size_t i, j;
public:
  struct value_type {
    std::size_t i, j;
  };

  tridiag_iterator( std::size_t i_, std::size_t j_ ) : i( i_ ), j( j_ ) { }

  value_type operator*() const {
    return { i, j };
  }

  bool operator!=( const tridiag_iterator &other ) {
    if( i >= other.i && j > other.j ) return 0;
    else return 1;
  }

  tridiag_iterator &operator++() {
    if( j == i ) j = i + 1;
    else if( j == i - 1 ) j = i;
    else if( j == i + 1 ) {
      i = i + 1;
      j = i - 1;
    }
  }
};

and I added two members to OneBodyMatrixElements:

  static tridiag_iterator begin() {
    return tridiag_iterator( 0, 0 );
  }

  static tridiag_iterator end() {
    return tridiag_iterator( N - 1, N - 1 );
  }

with N defined elsewhere. Usage:

  for( auto it = OneBodyMatrixElements::begin(), end = OneBodyMatrixElements::end(); it != end; ++it ) {
    OneBodyMatrixElements::val( (*it).i, (*it).j ) * func( (*it).i, (*it).j );
  }

Any comments are still welcome.

Answer 1

You could have OneBodyMatrixElements provide an iterator over the non-zero elements. Something like this:

struct Iterator
{
  struct value_type
  {
    std::size_t i, j;
    OneBodyMatrixElements::value_type val;
  };

  Iterator() : i(7), j(7) {}
  Iterator(std::size_t i, std::size_t j) : i(i), j(j) {}

  value_type operator* () const
  {
    return {i, j, OneBodyMatrixElements::val(i, j)};
  }

  Iterator& operator++ ()
  {
    if (j == i || j == i - 1) j = i + 1;
    else if (j == i + 1) { i = i +1; j = i - 1; }
    else :::
  }

private:
  std::size_t i, j;
};


struct OneBodyMatrixElements
{
  static Iterator begin()
  {
    return Iterator(0, 0);
  }

  static Iterator end()
  {
    return Iterator();
  }
};

// Usage:

for (auto it = OneBodyMatrixElements::begin(), end = OneBodyMatrixElements::end(); it != end; ++it)
{
  it->val * func(it->i, it->j);
}

The code above is intended to illustrate the idea, not provide a 100% ready solution. Other operator definitions, typedefs for the iterator, etc. have to be added.

Answer 2

My first approach would be to simply let the compiler figure out which ones are zero. Since your generator function is relatively (see below) simple and doesn't depend on any state outside itself, it should be feasible to analyze it far enough to inline it and to elide the unnecessary iterations.

My second approach would be to provide a cursor/iterator that only returns the non-zero element positions, possibly alongside with their values. In any case, profiling should be done first!

Anyhow, the compile might need some help:

• You static_cast the result ret, which already has that type, drop that. Instead, you might have to cast the assignments in the if-clauses.
• If all possible types can be used in a boolean context, use combined assignment and check to avoid an extra call: if (auto v = val(i, j)) ....
• Use else to make clear that the first and only the first if-clause that matches determines the result. Alternatively, use return instead of assigning to a temporary to make that clear.
• Another thing to think about your goal not to store the matrix. If you just spelled it out or computed it once, you would only have a few indexing operations that can usually be delivered from cache. Maybe the matrix in memory is even smaller than the code to compute it even!