Source

cexp / cexp / ext / src / constinloop.c

Full commit
typedef struct constinloop_{
  int num_s, num_n;
  double **x, **a, *b, e;
} ConstInLoop;

/* Dynamics:
     x(s+1) = (1 - e) * x(s) + e * (a x(t) + b)
 */

inline double dot(double * restrict a, double * restrict b, int num)
{
  int i;
  double c = 0;
  for (i = 0; i < num; ++i){
    c += a[i] * b[i];
  }
  return c;
}

inline void copy(double * restrict dest, double * restrict trg, int num)
{
  int i;
  for (i = 0; i < num; ++i){
    dest[i] = trg[i];
  }
}

/*[[[cog
import cog

template = """
inline void %(method)s(
    double ** restrict x, double ** restrict a,
    double * restrict%(const)s b, double%(const)s e,
    int%(const)s num_s, int%(const)s num_n)
{
  int s, n;
  for (s = 1; s < num_s; ++s){
    for (n = 0; n < num_n; ++n){
      x[s][n] += (1 - e) * x[s-1][n] + e * (dot(a[n], x[s-1], num_n) + b[n]);
    }
  }
}

int ConstInLoop_run_%(method)s(ConstInLoop *self, int repeat)
{
  int i;
  for(i = 0; i < repeat; ++i){
    %(method)s(self->x, self->a, self->b, self->e, self->num_s, self->num_n);
  }
  copy(self->x[0], self->x[self->num_s - 1], self->num_n);
  return 0;
}
"""

for (method, const) in [('with_const', ' const'), ('without_const', '')]:
    cog.out(template % dict(method=method, const=const))
]]]*/

inline void with_const(
    double ** restrict x, double ** restrict a,
    double * restrict const b, double const e,
    int const num_s, int const num_n)
{
  int s, n;
  for (s = 1; s < num_s; ++s){
    for (n = 0; n < num_n; ++n){
      x[s][n] += (1 - e) * x[s-1][n] + e * (dot(a[n], x[s-1], num_n) + b[n]);
    }
  }
}

int ConstInLoop_run_with_const(ConstInLoop *self, int repeat)
{
  int i;
  for(i = 0; i < repeat; ++i){
    with_const(self->x, self->a, self->b, self->e, self->num_s, self->num_n);
  }
  copy(self->x[0], self->x[self->num_s - 1], self->num_n);
  return 0;
}

inline void without_const(
    double ** restrict x, double ** restrict a,
    double * restrict b, double e,
    int num_s, int num_n)
{
  int s, n;
  for (s = 1; s < num_s; ++s){
    for (n = 0; n < num_n; ++n){
      x[s][n] += (1 - e) * x[s-1][n] + e * (dot(a[n], x[s-1], num_n) + b[n]);
    }
  }
}

int ConstInLoop_run_without_const(ConstInLoop *self, int repeat)
{
  int i;
  for(i = 0; i < repeat; ++i){
    without_const(self->x, self->a, self->b, self->e, self->num_s, self->num_n);
  }
  copy(self->x[0], self->x[self->num_s - 1], self->num_n);
  return 0;
}
/*[[[end]]]*/