Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Faulty Data Dependence Analysis #2194

Closed
AndrewCheng827 opened this issue Aug 31, 2023 · 5 comments · Fixed by #2212
Closed

Faulty Data Dependence Analysis #2194

AndrewCheng827 opened this issue Aug 31, 2023 · 5 comments · Fixed by #2212
Labels

Comments

@AndrewCheng827
Copy link
Contributor

AndrewCheng827 commented Aug 31, 2023

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
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Contributor

georgebisbas commented Sep 5, 2023

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
Copy link
Contributor

        (['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

@georgebisbas georgebisbas linked a pull request Sep 20, 2023 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
3 participants