Compile time array indexing using expression templates -- constexpr?

Question

I want to generate a function like this:

double apply_stencil(const double *u, const int i, const width,
                     const int *offset, const double *weight)
{
  double t=0;
  for(int j=0; j<width; j++)
    t += u[i+offset[j]] * weight[j];
  return t;
}

But I want to make sure that the width, the offsets, and possibly even the weights are compile/time constants.

That is possible to achieve by defining a type:

template <typename S>
double apply_stencil(const double *u, const int i)
{
  double t=0;
  for(int j=0; j<S::width; j++)
    t += u[i+S::offset[j]] * S::weight[j];
  return t;
}

// then:
struct Stencil {
  const static double weight[];
  const static int offset[];
  const static unsigned int width = 3;
};

const double Stencil::weight[] = {1.0, 2.0, 1.0};
const int Stencil::offset[]    = {-1, 0, 1};

However, this is not very pretty. I want a user to be able to specify the Stencil in their application code, and then call my apply_stencil function from my header file (this is really a simplification of something much more complicated).

Ideally, I would like to have the thing specified using expression templates, like so:

const Symbol u;
const Stencil<3> s (1*u[-1] + 2*u[0] + 1*u[1]);

Which uses this infrastructure:

struct Accessor {
  int offset;
};

struct Symbol
{
  Accessor operator[](int i) const {
    return Accessor{i};
  }
};

struct WeightedTerm
{
  double weight;
  int offset;
  WeightedTerm()
    : weight(1), offset(0) {}

  WeightedTerm(double a, Accessor u)
    : weight(a), offset(u.offset) {}
};

WeightedTerm operator*(double a, Accessor u) {
  return WeightedTerm(a,u);
}

template <int n>
struct Sum
{
  WeightedTerm terms[n];

  Sum(WeightedTerm x, WeightedTerm y) {
    terms[0] = x;
    terms[1] = y;
  }

  Sum(Sum<n-1> x, WeightedTerm y) {
    for(int i = 0; i < n-1; ++i)
      terms[i] = x.terms[i];

    terms[n-1] = y;
  }
};

Sum<2> operator+(WeightedTerm x, WeightedTerm y) {
  return Sum<2>(x,y);
}

template <int n>
Sum<n+1> operator+(Sum<n> x, WeightedTerm y) {
  return Sum<n+1>(x,y);
}

template <int width>
struct Stencil
{
  double weights[width];
  int offsets[width];
  Stencil(const Sum<width> s) {

    for(int j = 0; j < width; ++j) {
      weights[j] = s.terms[j].weight;
      offsets[j] = s.terms[j].offset;
    }
  };
};

This looks nice, but now, the arrays are not necessarily compile time known. If I write it like I do above with literals in the expression, I have verified that the compiler can make the correct optimizations. But I want to make find a way to guarantee that they are always compile time constants.

I presume I could encode the offset as a template parameter in Accessor and WeightedTerm, but I can't see how I could do that and keep the desired expression syntax since operator() takes the offset as a regular argument.

So, the question is if there is a way to achieve this? I have a feeling that constexpr could be of use here which I am a bit unfamiliar with.

Yakk - Adam Nevraumont · Accepted Answer · 2016-12-30 16:22:56Z

1

Use a std::integral_constant alias. If you want ease of use, use user defined literals.

constexpr int int_from_chars(){return 0;}

constexpr int k_pow(unsigned int x, unsigned int b=10){
  return (x==0)?1:(k_pow(x-1,b)*b);
}

template<class...Chars>
constexpr int int_from_chars(char x, Chars...chars){
  return k_pow(sizeof...(Chars))*(x-'0') + int_from_chars(chars...);
}
template<char...Cs>
constexpr auto operator""_kint(){
  return std::integral_constant<int, int_from_chars(Cs...)>{};
}

or somesuch. Needs polish.

Now 7_kint is a compile time type encoding 7.

You can take the integral_constant<int, K> as a function argument, where K is a template non type argument.

Extending to octal/hex/binary left as an exercise.

template<int width>
double apply_stencil(const double *u, const int i, std::integral_constant<int,width>,
                 const int *offset, const double *weight)
{
  double t=0;
  for(int j=0; j<width; j++)
    t += u[i+offset[j]] * weight[j];
  return t;
}

Called via

apply_stencil( ptr, i, 7_kint, poffsets, pweights);

Similar (but more complex) techniques can be used to pass packs of offsets. Weights are a mess, as there is little compile time double support.

edited Dec 30, 2016 at 16:22

answered Dec 30, 2016 at 15:38

Yakk - Adam Nevraumont

280k31 gold badges367 silver badges585 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

kalj Over a year ago

Interesting. I assume this would be used then like so: const Stencil<3> s (1*u[m1] + 2*u[z] + 1*u[p1]);. Correct? Also, "user defined literals" -- what do you refer to by that?

Yakk - Adam Nevraumont Over a year ago

@kalj please google the term before askimg about it. You woukd pass the _intk to apply_stencil as arguments, but their value would be guaranteed compile time.

Collectives™ on Stack Overflow

Compile time array indexing using expression templates -- constexpr?

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related