ConditionalDimension: Devito produces a code referencing a dimension before its declaration

[maksgraczyk]

When ConditionalDimension is used with the intended inner-most dimension as its parent and more than one dimension in the condition, Devito sometimes produces a code referencing one of the dimensions before its declaration. This leads to a compilation error.

MFE:

from devito import *
from sympy import And, Or

t1, t2, t3, t4, t5, t6, t7, t8, t9, t10 = \
    dimensions('t1 t2 t3 t4 t5 t6 t7 t8 t9 t10')

f0 = Function(name='f0', grid=Grid(shape=(2, 2, 4, 4),
                                   dimensions=(t1, t2, t3, t4)))
f1 = Function(name='f1', grid=Grid(shape=(2, 2, 3, 3),
                                   dimensions=(t5, t6, t7, t8)))
f2 = Function(name='f2', grid=f1.grid)

cd = ConditionalDimension(name='cd', parent=t10,
                          condition=Or(Gt(f0[t5, t6, t7 + t9,
                                             t8 + t10],
                                          f1[t5, t6, t7, t8]),
                                       And(~Ne(f0[t5, t6, t7 + t9,
                                                  t8 + t10],
                                               f1[t5, t6, t7, t8]),
                                           Lt(2 * t9 + t10,
                                              f2[t5, t6, t7, t8]))))

op = Operator([Eq(f1[t5, t6, t7, t8], f0[t5, t6, t7 + t9, t8 + t10],
                  implicit_dims=cd),
               Eq(f2[t5, t6, t7, t8], 2 * t9 + t10, implicit_dims=cd)])

op.apply(t9_M=1, t10_M=1)

Traceback:

/tmp/devito-jitcache-uid1000/54ae8b2b7dd776276090529a8f3645f1098562ce.c: In function ‘Kernel’:
/tmp/devito-jitcache-uid1000/54ae8b2b7dd776276090529a8f3645f1098562ce.c:48:43: error: ‘t8’ undeclared (first use in this function); did you mean ‘t9’?
   48 |             if (f0[t5][t6][t7 + t9][t10 + t8] == f1[t5][t6][t7][t8] && t10 + 2*t9 < f2[t5][t6][t7][t8] || f0[t5][t6][t7 + t9][t10 + t8] > f1[t5][t6][t7][t8])
      |                                           ^~
      |                                           t9
/tmp/devito-jitcache-uid1000/54ae8b2b7dd776276090529a8f3645f1098562ce.c:48:43: note: each undeclared identifier is reported only once for each function it appears in
FAILED compiler invocation:gcc -O3 -g -fPIC -Wall -std=c99 -march=native -Wno-unused-result -Wno-unused-variable -Wno-unused-but-set-variable -ffast-math -shared -fopenmp /tmp/devito-jitcache-uid1000/54ae8b2b7dd776276090529a8f3645f1098562ce.c -lm -o /tmp/devito-jitcache-uid1000/54ae8b2b7dd776276090529a8f3645f1098562ce.so
Traceback (most recent call last):
  File "mfe.py", line 27, in <module>
    op.apply(t9_M=1, t10_M=1)
  File "/home/maksymilian/Desktop/UROP/devito/devito/operator/operator.py", line 674, in apply
    cfunction = self.cfunction
  File "/home/maksymilian/Desktop/UROP/devito/devito/operator/operator.py", line 589, in cfunction
    self._jit_compile()
  File "/home/maksymilian/Desktop/UROP/devito/devito/operator/operator.py", line 574, in _jit_compile
    recompiled, src_file = self._compiler.jit_compile(self._soname,
  File "/home/maksymilian/Desktop/UROP/devito/devito/compiler.py", line 316, in jit_compile
    _, _, _, recompiled = compile_from_string(self, target, code, src_file,
  File "/home/maksymilian/.conda/envs/devito/lib/python3.8/site-packages/codepy/jit.py", line 433, in compile_from_string
    toolchain.build_extension(ext_file, source_paths, debug=debug)
  File "/home/maksymilian/.conda/envs/devito/lib/python3.8/site-packages/codepy/toolchain.py", line 210, in build_extension
    raise CompileError("module compilation failed")
codepy.CompileError: module compilation failed

Generated code:

#define _POSIX_C_SOURCE 200809L
#include "stdlib.h"
#include "math.h"
#include "sys/time.h"
#include "xmmintrin.h"
#include "pmmintrin.h"

struct dataobj
{
  void *restrict data;
  int * size;
  int * npsize;
  int * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
} ;

struct profiler
{
  double section0;
} ;


int Kernel(struct dataobj *restrict f0_vec, struct dataobj *restrict f1_vec, struct dataobj *restrict f2_vec, const int t10_M, const int t10_m, const int t5_M, const int t5_m, const int t6_M, const int t6_m, const int t7_M, const int t7_m, const int t8_M, const int t8_m, const int t9_M, const int t9_m, struct profiler * timers)
{
  float (*restrict f0)[f0_vec->size[1]][f0_vec->size[2]][f0_vec->size[3]] __attribute__ ((aligned (64))) = (float (*)[f0_vec->size[1]][f0_vec->size[2]][f0_vec->size[3]]) f0_vec->data;
  float (*restrict f1)[f1_vec->size[1]][f1_vec->size[2]][f1_vec->size[3]] __attribute__ ((aligned (64))) = (float (*)[f1_vec->size[1]][f1_vec->size[2]][f1_vec->size[3]]) f1_vec->data;
  float (*restrict f2)[f2_vec->size[1]][f2_vec->size[2]][f2_vec->size[3]] __attribute__ ((aligned (64))) = (float (*)[f2_vec->size[1]][f2_vec->size[2]][f2_vec->size[3]]) f2_vec->data;

  /* Flush denormal numbers to zero in hardware */
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
  struct timeval start_section0, end_section0;
  gettimeofday(&start_section0, NULL);
  /* Begin section0 */
  for (int t5 = t5_m; t5 <= t5_M; t5 += 1)
  {
    for (int t6 = t6_m; t6 <= t6_M; t6 += 1)
    {
      for (int t7 = t7_m; t7 <= t7_M; t7 += 1)
      {
        for (int t9 = t9_m; t9 <= t9_M; t9 += 1)
        {
          #pragma omp simd aligned(f0,f1,f2:32)
          for (int t10 = t10_m; t10 <= t10_M; t10 += 1)
          {
            if (f0[t5][t6][t7 + t9][t10 + t8] == f1[t5][t6][t7][t8] && t10 + 2*t9 < f2[t5][t6][t7][t8] || f0[t5][t6][t7 + t9][t10 + t8] > f1[t5][t6][t7][t8])
            {
              for (int t8 = t8_m; t8 <= t8_M; t8 += 1)
              {
                f1[t5][t6][t7][t8] = f0[t5][t6][t7 + t9][t10 + t8];
                f2[t5][t6][t7][t8] = t10 + 2*t9;
              }
            }
          }
        }
      }
    }
  }
  /* End section0 */
  gettimeofday(&end_section0, NULL);
  timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
  return 0;
}

[mloubout]

None of f0, f1 and f2 have t9 or t10 as a dimension so not sure it's supposed to work since the compiler may not be able to understand it will be used as one

[mloubout]

Also I'd like to point out that the condition here is non-unique you can have multuple pairs of t8, t10 such that t8 + t10 is the same value (in a symbolic interpretation way which is what the compiler does, not look at actual bounds of each). So a single parent dimension can't figure that out since it depends on two dimensions put together.

So that may be a bug but that something is may not be easily automatable since it depends on the actual values of for example t8_m, t8_M, t10_m, t10M to know where best to put the condition.

In the mean time, isn't there a way avoid these two extra dimension and use Constant and/or t8.spacing/t7.spacing?

[FabioLuporini]

this was fixed by the recent #2007 , hence closing

I've already added the MFE to an upcoming PR

[georgebisbas]

Maks' mfe not working for me?

    f0 = Function(name='f0', grid=Grid(shape=(2, 2, 4, 4),
  File "/home/gb4018/workspace/devito/devito/types/grid.py", line 157, in __init__
    raise ValueError("Cannot create Grid with Dimension `%s` "
ValueError: Cannot create Grid with Dimension `t1` since it's not a SpaceDimension

[FabioLuporini]

Have u actually read the error message?

[georgebisbas]

Yes all good, I was only saying that the above specific MFE is not a copy paste in the tests