diku-dk / futhark

:boom::computer::boom: A data-parallel functional programming language
http://futhark-lang.org
ISC License
2.37k stars 164 forks source link

pyopencl --library : Internal compiler error #763

Closed LaurieGustavo closed 5 years ago

LaurieGustavo commented 5 years ago

Hello,

i have an internal compiler error for my program when trying to compile to a Python module. Would anybody have an idea about a fix ? thank you =)

it's a program for the Lattice-Baltzmann algorithm. Basically there are 4 matrices that need each other to be updated at each iteration (u: [nx][ny](f64, f64), rho: [nx][ny]f64, f_in: [nx][ny][9], f_out: [nx][ny][9]f64 ). The problem seems to be that it can't determine the size of f_in or f_out

futhark --version Futhark 0.11.0 git: cb04266 (Mon May 6 14:24:26 2019 +0200)

futhark pyopencl --library Lattice_Boltzmann.fut
Internal compiler error.
Please report this at https://github.com/diku-dk/futhark/issues.
Type error after pass 'simplify':
In function main
When checking function body
In expression of statement {[num_elems_20529][num_elems_20531][9i32]f64 res_20539}
Type error:
In function main
When checking function body
In expression of statement {[num_elems_20529][num_elems_20531][9i32]f64 res_20539}
When checking certificates
In subexp index_certs_20911
Use of unknown variable index_certs_20911.

Lattice_Baltzmann.fut

-- Program for Lattice-Boltzmann in Futhark

import "vector"

------------------------------------
----- FONCTIONS INITIALISATION -----
------------------------------------

-- Initilization of a matrix nx * ny with value x
-- parametric polymorphism
let init_matrix 't (nx: i32)(ny: i32)(x: t): [nx][ny]t =
    map( \(_):[]t -> 
        map( \(_):t ->
            x
        ) (0..<ny)
    ) (0..<nx)

-- Initilization of an array nx * 1 with value x
-- parametric polymorphism
let init_array 't (nx: i32)(x: t): [nx]t =
    map( \(_):t ->
        x
    ) (0..<nx)

------------------------------------
------------ CONSTANTES ------------
------------------------------------

let c: [9]vector        = [ ( 0,  0),
                            (-1,  1),
                            (-1,  0),
                            (-1, -1),
                            ( 0, -1),
                            ( 1, -1),
                            ( 1,  0),
                            ( 1,  1),
                            ( 0,  1) ]

let opposite: [9]i32    = [0, 5, 6, 7, 8, 1, 2, 3, 4]

let w: [9]f64           = [4/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9]

------------------------------------
------- FONCTIONS DE CALCULS -------
------------------------------------

-- calculates rho at some point
let rho_calcul(f_in: []f64, mask: u32): f64 =
    if mask > 0 then 0 else reduce (+) 0 f_in

-- updates rho at every point
let rho_update [nx][ny] (f_in: [nx][ny][9]f64, mask: [nx][ny]u32): [nx][ny]f64 =
    map2(\fx mx ->
        map2(\fxy mxy ->
            rho_calcul(fxy, mxy)
        ) fx mx
    ) f_in mask

-- calculates u at point (x,y)
let u_calcul(f_in: []f64, rho: f64, mask: u32): vector =
    if mask > 0 then
        (0,0)
    else
        let res = reduce add (0,0) (map2 mult f_in c)
        in div(rho,res)

-- updates u at every point
let u_update [nx][ny](f_in: [nx][ny][9]f64, rho: [nx][ny]f64, mask: [nx][ny]u32): [nx][ny]vector =
    map3 (\fx rx mx ->
        map3 ( \fxy rxy mxy ->
            u_calcul(fxy, rxy, mxy)
        ) fx rx mx
    ) f_in rho mask

-- j(x,y) = u[x][y] + Tau*g
let j(u: vector, tau: f64, g: vector): vector =
    let a = mult tau g
    in add a u

-- f_eq
let f_eq(i: i32, rho: f64, u: vector, g: vector, tau: f64): f64 =
    let j_val   = j(u,tau,g)
    let dot_val = dotprod(c[i],j_val)
    in w[i]*rho*(1 + 3*dot_val + (9f64/2)*dot_val*dot_val - (3f64/2)*square(j_val))

-- Initilization of f_in
let init_f_in [nx][ny] (rho: [nx][ny]f64, u: [nx][ny]vector, g: vector, tau: f64): [nx][ny][9]f64 =
    map( \(x) -> 
        map( \(y) ->
            map( \(i) ->
                f_eq(i, rho[x,y], u[x,y], g, tau)
            ) (0..<9)
        ) (0..<ny)
    ) (0..<nx)

-- calculates f_out at point (x,y,i)
let f_out_calcul(f_in: []f64, i: i32, rho: f64, u: vector, g: vector, tau: f64, mask: u32): f64 =   

    if mask > 0 then
        f_in[opposite[i]]
    else
        f_in[i] - (f_in[i] - f_eq(i,rho,u,g,tau))/tau 

-- updates f_out at every point
let f_out_update [nx][ny] (f_in: [nx][ny][9]f64, rho: [nx][ny]f64, u: [nx][ny]vector, tau: f64, g: vector, mask: [nx][ny]u32): [nx][ny][9]f64 =
    map( \(x) -> 
        map( \(y) ->
            map( \(i) ->
                f_out_calcul(f_in[x,y], i, rho[x,y], u[x,y], g, tau, mask[x,y])
            ) (0..<9)
        ) (0..<ny)
    ) (0..<nx)

-- collision phase
let collision [nx][ny] (f_in: [nx][ny][9]f64, g: vector, tau: f64, mask: [nx][ny]u32): [nx][ny][9]f64 =
    let rho     = rho_update(f_in, mask)
    let u       = u_update(f_in, rho, mask)
    in f_out_update(f_in, rho, u, tau, g, mask)

-- propagation phase
let propagation [nx][ny] (f_out: [nx][ny][9]f64): [nx][ny][9]f64 =
    let f_in = init_matrix(nx)(ny)(init_array(9)(0f64))
    in loop (f_in) for x < nx do
        loop (f_in) for y < ny do
            loop (f_in) for i < 9 do 
                f_in with [(x+(i32.f64 c[i].1)+nx)%nx,(y+(i32.f64 c[i].2)+ny)%ny,i] = f_out[x,y,i]

let norm_u [nx][ny] (u: [nx][ny]vector): [nx][ny]f64 =
    map (\(u_x) ->
        map (\(u_x_y) ->
            norm(u_x_y)
        ) u_x
    ) u

-- Lattice-Boltzmann
let main (nx: i32)(ny: i32)(g_x: f64)(g_y: f64)(tau: f64)(it: i32): [nx][ny]f64 =

    -- Constantes
    let g: vector   = (g_x, g_y)
    let mask        = init_matrix(nx)(ny)(0u32)

    -- Initialize u
    let u           = init_matrix(nx)(ny)((0f64, 0f64))
    -- Initialize rho
    let rho         = init_matrix(nx)(ny)(1f64)

    -- Initialize f_in & f_out
    let f_in        = init_f_in(rho, u, g, tau)
    let f_out       = f_out_update(f_in, rho, u, tau, g, mask)
    let f_in        = propagation(f_out)

    let (f_in) = loop (f_in) for _i < it do
        let f_out = collision(f_in, g, tau, mask)
        in (propagation(f_out))

    let rho         = rho_update(f_in, mask)
    let u           = u_update(f_in, rho, mask)

    in norm_u(u)

vector.fut


-- VECTOR file --

import "futlib/math"

-- type vector
type vector = (f64,f64)

------------ FUNCTIONS ------------

-- addition of two vectors
let add(v1: vector)(v2: vector): vector =
    let(a,b) = v1
    let(c,d) = v2
    in (a+c, b+d)

-- substract a vector to another
let sub(v1: vector, v2: vector): vector =
    let(a,b) = v1
    let(c,d) = v2
    in (a-c, b-d)

-- multiply a vector by a float
let mult(f: f64)(v: vector): vector =
    let (a,b) = v
    in (f*a, f*b)

-- divide a vector by a float
let div(f: f64, v: vector): vector =
    let(a,b) = v
    in (a/f, b/f)

-- dotprod of two vectors
let dotprod(v1: vector, v2: vector): f64 =
    let(a,b) = v1
    let(c,d) = v2
    in a*c + b*d

-- square of a vector
let square(v: vector): f64 =
    dotprod(v,v)

-- norm of a vector
let norm(v: vector): f64 =
    f64.sqrt(square(v))
athas commented 5 years ago

Thanks for the report; I'll take a look. I don't expect this one will be too difficult to fix. As a workaround, put unsafe the beginning of the main function (with no bounds checking, there will be no safety certificates that the compiler can screw up).

athas commented 5 years ago

Smallest reproducible case:

type vector = (f64,f64)

let add(v1: vector)(v2: vector): vector =
    let(a,b) = v1
    let(c,d) = v2
    in (a+c, b+d)

let mult(f: f64)(v: vector): vector =
    let (a,b) = v
    in (f*a, f*b)

let dotprod(v1: vector, v2: vector): f64 =
    let(a,b) = v1
    let(c,d) = v2
    in a*c + b*d

let square(v: vector): f64 =
    dotprod(v,v)

let init_matrix 't (nx: i32)(ny: i32)(x: t): [nx][ny]t =
    map( \(_):[]t ->
        map( \(_):t ->
            x
        ) (0..<ny)
    ) (0..<nx)

let j(u: vector, tau: f64, g: vector): vector =
    let a = mult tau g
    in add a u

let f_eq(rho: f64, u: vector, g: vector, tau: f64): f64 =
    let j_val   = j(u,tau,g)
    in rho*(1 + 3 + (9f64/2) - (3f64/2)*square(j_val))

let init_f_in [nx][ny] (rho: [nx][ny]f64, u: [nx][ny]vector, g: vector, tau: f64): [nx][ny][9]f64 =
    map( \(x) ->
        map( \(y) ->
                       replicate 9 (f_eq(rho[x,y], u[x,y], g, tau))
        ) (0..<ny)
    ) (0..<nx)

let main (nx: i32)(ny: i32)(g_x: f64)(g_y: f64)(tau: f64) =
    let g: vector   = (g_x, g_y)
    let u           = init_matrix(nx)(ny)((0f64, 0f64))
    let rho         = init_matrix(nx)(ny)(1f64)
    let f_in        = init_f_in(rho, u, g, tau)
    in f_in
LaurieGustavo commented 5 years ago

it seems the error disappears when using map2 instead of map , for example in function init_f_in in the smallest reproducible case :

let init_f_in [nx][ny] (rho: [nx][ny]f64, u: [nx][ny]vector, g: vector, tau: f64): [nx][ny][9]f64 =
    map2(\rx ux ->
        map2(\rxy uxy ->
            replicate 9 (f_eq(rxy, uxy, g, tau))
        ) rx ux
    ) rho u

So know i can compile my main code to a Python module =)

athas commented 5 years ago

Have you tried the newest release of the compiler? I am pretty sure I fixed this even for your original formulation.