Faulty Data Dependence Analysis

[AndrewCheng827]

Minimal Failing Example:

from devito import *

def solver():  
    grid = Grid(shape=(10,))    
    u = TimeFunction(name='u',grid=grid)

    x = grid.dimensions[0]
    h_x = x.spacing

    eq_1 = Eq(u[0, x], 1)
    eq_3 = Eq(u[1,x], u[0, x + h_x] + u[0, x - h_x] - 2*u[0,x])
    op = Operator([eq_1, eq_3], opt='noop')

    print("C CODE: ")
    print(op.ccode)

if __name__ == '__main__':
    solver()

Incorrect output (C Code):

int Kernel(struct dataobj *restrict u_vec, const int x_M, const int x_m, struct profiler * timers)
{
  float (*restrict u)[u_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[u_vec->size[1]]) u_vec->data;

  /* Begin section0 */
  START_TIMER(section0)
  for (int x = x_m; x <= x_M; x += 1)
  {
    u[0][x + 1] = 1;

    u[1][x + 1] = u[0][x] - 2*u[0][x + 1] + u[0][x + 2];
  }
  STOP_TIMER(section0,timers)
  /* End section0 */

  return 0;
}

Expected Output:
[[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[-1. 0. 0. 0. 0. 0. 0. 0. 0. -1.]]

Actual Output:
[[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[-2. -1. -1. -1. -1. -1. -1. -1. -1. -1.]]

Issue:
The second line in the loop u[1][x+1] = u[0][x] - 2*u[0][x+1] + u[0][x+2] has a dependency on [u][0][x+2]. However that value is not initialized yet at the time the above equation is computed.

Bisbas edit:

New MFE added: #2194 (comment)

[georgebisbas]

Modelling the same thing as well (not mfe, but better way to model):

import numpy as np
from devito import Grid, TimeFunction, Eq, Operator


def solver():
    grid = Grid(shape=(10, ))
    u = TimeFunction(name='u', grid=grid)

    x = grid.dimensions[0]
    h_x = x.spacing
    t = grid.stepping_dim
    
    u.data[:, :] = 1
    eq0 = Eq(u[t, x], u[t - 1, x + h_x] - 2*u[t - 1, x] + u[t - 1, x - h_x])

    op = Operator([eq0], opt='noop')
    op.apply(time_M=1)
    # print("C CODE: ")
    # print(op.ccode)

    expected = np.array([[ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
                         [-1., 0., 0., 0., 0., 0., 0., 0., 0., -1.]])

    assert(np.all(u.data[:] == expected[:]))


if __name__ == '__main__':
    solver()

[georgebisbas]

or

import numpy as np
from devito import Grid, TimeFunction, Eq, Operator


def solver():
    grid = Grid(shape=(10, ))
    u = TimeFunction(name='u', grid=grid)

    x = grid.dimensions[0]
    h_x = x.spacing
    t = grid.stepping_dim
    
    eq_1 = Eq(u[t - 1, x], 1)
    eq_3 = Eq(u[t, x], u[t - 1, x + h_x] - 2*u[t - 1, x] + u[t - 1, x - h_x])

    op = Operator([eq_1, eq_3], opt='noop')
    op.apply(time_M=1)
    # print("C CODE: ")
    # print(op.ccode)

    expected = np.array([[ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
                         [-1., 0., 0., 0., 0., 0., 0., 0., 0., -1.]])

    assert(np.all(u.data[:] == expected[:]))


if __name__ == '__main__':
    solver()

[georgebisbas]

Initial with test:

import numpy as np
from devito import Grid, TimeFunction, Eq, Operator


def solver():
    grid = Grid(shape=(10, ))
    u = TimeFunction(name='u', grid=grid)

    x = grid.dimensions[0]
    h_x = x.spacing

    eq_1 = Eq(u[0, x], 1)
    eq_3 = Eq(u[1, x], u[0, x + h_x] - u[0, x])

    op = Operator([eq_1, eq_3])
    op.apply()
    print("C CODE: ")
    print(op.ccode)

    expected = np.array([[ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.], 
                         [-1., 0., 0., 0., 0., 0., 0., 0., 0., -1.]])

    assert(np.all(u.data[:] == expected[:]))


if __name__ == '__main__':
    solver()

[georgebisbas]

New MFE:

import numpy as np
from devito import Grid, TimeFunction, Eq, Operator, Dimension, Function

grid = Grid(shape=(2, 2))

u0 = Function(name='u0', grid=grid)
x, y = grid.dimensions
t = grid.stepping_dim

eq0 = Eq(u0[0, y], 1)
eq1 = Eq(u0[0+1, y], u0[0, y + 1])

op = Operator([eq0, eq1])
print(op.ccode)
op.apply()

expected = np.array([[ 1., 1.],
                     [ 1., 1.]])
import pdb;pdb.set_trace()
assert(np.all(u0.data[:] == expected[:]))

[georgebisbas]

        (['Eq(u[0, y], 1)', 'Eq(u[1, 1], u[0, y + 1])'],
            np.array([[1., 1.], [0., 0.]]),
            ['y', 'y'], 'y,y'),
        (['Eq(u[0, 1], 1)', 'Eq(u[x, y], u[0, y])'],
            np.array([[0., 1.], [0., 1.]]),
            ['xy'], 'x,y')
    ])
    def test_2194_v2(self, eqns, expected, exp_trees, exp_iters):
        grid = Grid(shape=(2, 2))
        u = Function(name='u', grid=grid)
        x, y = grid.dimensions

        for i, e in enumerate(list(eqns)):
            eqns[i] = eval(e)

        op = Operator(eqns)
        assert_structure(op, exp_trees, exp_iters)

        op.apply()
        assert(np.all(u.data[:] == expected[:]))

This test breaks with openmp, should look better tmr and maybe open new issue