Work in progress: Initial Haskell version of MUDA

by syoyo

[En]

Finally, I’ve finished initial Haskell reimplementation of MUDA language[1,2, 3] compiler.

http://lucille.svn.sourceforge.net/viewvc/lucille/angelina/haskellmuda/

Here is an example.


// input.mu
// Input MUDA code
vec add_func(vec a, vec b)
{
        return a + b;
}

From this input, MUDA compiler(with x86 SSE CodeGen) generates the following SIMD code.


$ mudah input.mu

//
// The following code is generated by MUDA compiler
//
#include 
#ifdef __64bit__
    #error "64bit env is not yet supported"
#else
#endif

#ifdef __GNUC__
static inline void *muda_aligned_addr16(void *addr) {
    return (void *)((((unsigned int)addr) + 15UL) & ~(15UL));
}
#endif

#ifdef __GNUC__
    #define MUDA_ATTRIB_ALIGN __attribute__((aligned(16)))
    #define MUDA_DECL_ALIGN
#elif defined(_MSC_VER)
    #define MUDA_ATTRIB_ALIGN
    #define MUDA_DECL_ALIGN __declspec(align(16))
#else
    #error "Sorry, MUDA doesn't support your compiler
#endif

static inline __m128 muda_sel_ps( __m128 a, __m128 b, __m128 mask )
{
    b = _mm_and_ps( b, mask );
    a = _mm_andnot_ps( mask, a);
    return _mm_or_ps( a, b );
}

__m128 add_func (__m128 a, __m128 b)
{
    const __m128 t_vec1  =  _mm_add_ps( a ,  b  ) ;
    return t_vec1 ;
}

And more, another example. SIMD ray-triangle intersection code written in MUDA.


// isect.mu
//
// SIMD cross product
static inline vec cross4(vec a, vec b, vec c, vec d)
{
        return ((a * b) - (c * d));
}

// SIMD dot product
static inline vec dot4(vec ax, vec ay, vec az, vec bx, vec by, vec bz)
{
        return (ax * bx + ay * by + az * bz);
}

//
// 1 ray - 4 triangles intersection.
// Assume that, for ray data, same value is copied over 4 elements.
//
vec
isect(vec rox, vec roy, vec roz,
      vec rdx, vec rdy, vec rdz,
      vec v0x, vec v0y, vec v0z,
      vec e1x, vec e1y, vec e1z,
      vec e2x, vec e2y, vec e2z,
      out vec outT, out vec outU, out vec outV, vec prevT)
{
        vec px, py, pz;

        // p = d x e2
        px = cross4(e2z, e2y, rdy, rdz);
        py = cross4(e2x, e2z, rdz, rdx);
        pz = cross4(e2y, e2x, rdx, rdy);

        vec sx, sy, sz;

        sx = rox - v0x;
        sy = roy - v0y;
        sz = roz - v0z;

        vec qx, qy, qz;

        vec vone  = vec(1.0);
        vec vzero = vec(0.0);

        vec a     = dot4(px, py, pz, e1x, e1y, e1z);

        // `slash dot' operator means approximate division if possible.
        vec rpa   = vone /. a;

        // q = s x e1
        qx = cross4(e1z, e1y, sy, sz);
        qy = cross4(e1x, e1z, sz, sx);
        qz = cross4(e1y, e1x, sx, sy);

        vec u, v, t;

        u = dot4(sx, sy, sz, px, py, pz) * rpa;
        v = dot4(rdx, rdy, rdz, qx, qy, qz) * rpa;
        t = dot4(e2x, e2y, e2z, qx, qy, qz) * rpa;

        vec eps;

        eps = vec(0.00001);

        vec mask0;

        mask0 = (((a * a) > eps) & ((u + v)  vzero) & (v > vzero));

        vec mask;

        //
        // final mask
        //
        mask = (mask0 & t) & (t < prevT);

        outT = sel(outT, t, mask);
        outU = sel(outU, t, mask);
        outV = sel(outV, t, mask);

        return mask;
}

MUDA compiler compiles this code into…


$ mudah isect.mu

...

__m128 cross4 (__m128 a, __m128 b, __m128 c, __m128 d)
{
    const __m128 t_vec1  =  _mm_mul_ps( a ,  b  ) ;
    const __m128 t_vec2  =  _mm_mul_ps( c ,  d  ) ;
    const __m128 t_vec3  =  _mm_sub_ps( t_vec1 ,  t_vec2  ) ;
    return t_vec3 ;
}

__m128 dot4 (__m128 ax, __m128 ay, __m128 az, __m128 bx, __m128 by, __m128 bz)
{
    const __m128 t_vec4  =  _mm_mul_ps( ax ,  bx  ) ;
    const __m128 t_vec5  =  _mm_mul_ps( ay ,  by  ) ;
    const __m128 t_vec6  =  _mm_add_ps( t_vec4 ,  t_vec5  ) ;
    const __m128 t_vec7  =  _mm_mul_ps( az ,  bz  ) ;
    const __m128 t_vec8  =  _mm_add_ps( t_vec6 ,  t_vec7  ) ;
    return t_vec8 ;
}

__m128 isect (__m128 rox, __m128 roy, __m128 roz,
__m128 rdx, __m128 rdy, __m128 rdz,
__m128 v0x, __m128 v0y, __m128 v0z,
__m128 e1x, __m128 e1y, __m128 e1z,
__m128 e2x, __m128 e2y, __m128 e2z,
__m128 * outT, __m128 * outU, __m128 * outV, __m128 prevT)
{
    __m128 px ;

    __m128 py ;

    __m128 pz ;

    const __m128 t_vec9 = cross4 (e2z, e2y, rdy, rdz) ;
    px = t_vec9 ;

    const __m128 t_vec10 = cross4 (e2x, e2z, rdz, rdx) ;
    py = t_vec10 ;

    const __m128 t_vec11 = cross4 (e2y, e2x, rdx, rdy) ;
    pz = t_vec11 ;

    __m128 sx ;

    __m128 sy ;

    __m128 sz ;

    const __m128 t_vec12  =  _mm_sub_ps( rox ,  v0x  ) ;
    sx = t_vec12 ;

    const __m128 t_vec13  =  _mm_sub_ps( roy ,  v0y  ) ;
    sy = t_vec13 ;

    const __m128 t_vec14  =  _mm_sub_ps( roz ,  v0z  ) ;
    sz = t_vec14 ;

    __m128 qx ;

    __m128 qy ;

    __m128 qz ;

    __m128 vone ;

    __m128 vzero ;

    __m128 a ;

    __m128 rpa ;

    const __m128 t_vec15 = cross4 (e1z, e1y, sy, sz) ;
    qx = t_vec15 ;

    const __m128 t_vec16 = cross4 (e1x, e1z, sz, sx) ;
    qy = t_vec16 ;

    const __m128 t_vec17 = cross4 (e1y, e1x, sx, sy) ;
    qz = t_vec17 ;

    __m128 u ;

    __m128 v ;

    __m128 t ;

    const __m128 t_vec18 = dot4 (sx, sy, sz, px, py, pz) ;
    const __m128 t_vec19  =  _mm_mul_ps( t_vec18 ,  rpa  ) ;
    u = t_vec19 ;

    const __m128 t_vec20 = dot4 (rdx, rdy, rdz, qx, qy, qz) ;
    const __m128 t_vec21  =  _mm_mul_ps( t_vec20 ,  rpa  ) ;
    v = t_vec21 ;

    const __m128 t_vec22 = dot4 (e2x, e2y, e2z, qx, qy, qz) ;
    const __m128 t_vec23  =  _mm_mul_ps( t_vec22 ,  rpa  ) ;
    t = t_vec23 ;

    __m128 eps ;

    const float t_float24 = 1.0e-5 ;
    const __m128 t_vec25 = _mm_set_ps1(  t_float24  ) ;
    const __m128 t_vec26 = t_vec25 ;
    eps = t_vec26 ;

    __m128 mask0 ;

    const __m128 t_vec27  =  _mm_mul_ps( a ,  a  ) ;
    const __m128 t_vec28  =  _mm_cmpgt_ps( t_vec27 ,  eps  ) ;
    const __m128 t_vec29  =  _mm_add_ps( u ,  v  ) ;
    const __m128 t_vec30  =  _mm_cmplt_ps( t_vec29 ,  vone  ) ;
    const __m128 t_vec31  =  _mm_and_ps( t_vec28 ,  t_vec30  ) ;
    const __m128 t_vec32  =  _mm_cmpgt_ps( u ,  vzero  ) ;
    const __m128 t_vec33  =  _mm_cmpgt_ps( v ,  vzero  ) ;
    const __m128 t_vec34  =  _mm_and_ps( t_vec32 ,  t_vec33  ) ;
    const __m128 t_vec35  =  _mm_and_ps( t_vec31 ,  t_vec34  ) ;
    mask0 = t_vec35 ;

    __m128 mask ;

    const __m128 t_vec36  =  _mm_and_ps( mask0 ,  t  ) ;
    const __m128 t_vec37  =  _mm_cmplt_ps( t ,  prevT  ) ;
    const __m128 t_vec38  =  _mm_and_ps( t_vec36 ,  t_vec37  ) ;
    mask = t_vec38 ;

    const __m128 t_vec39 = muda_sel_ps ((*outT), t, mask) ;
    (*outT) = t_vec39 ;

    const __m128 t_vec40 = muda_sel_ps ((*outU), t, mask) ;
    (*outU) = t_vec40 ;

    const __m128 t_vec41 = muda_sel_ps ((*outV), t, mask) ;
    (*outV) = t_vec41 ;

    return mask ;
}

TODOs

Here is TODOs I’m planning until the end of this year.

– Documentation!
– Support for struct
– Support for C/C++ integration(header file generation)
CodeGen for LLVM IR
– swizzling. e.g. vec dst = src.xxyy
– Support for integer and int x 4, double x 4 and double x 2 SIMD.
– Integration with AMD simulation utilities(For counting instructions)

[Ja]

とりあえず Haskell で MUDA コンパイラ書き直しました。
OCaml で初期実装したときの状態くらいにまではなっています。

ここまで時間がかかったのは、シンボルテーブルの管理でどうしても逐次処理 & 大域変数処理を伴うので、
それをどう Haskell にマッピングするかをうんうん考えていたからでした。

Hashtable, seq, unsafePerformIO とかいろいろ試したけど、結局のところ State モナドを使うことに。
Haskell 習得(主にモナド)も含めてですが、そんなんで一ヶ月半くらい費やしてしまった。

インド工科大(IIT)生なら一週間で C 言語を習得するとのことですから、
私ももっと精進しないとですね。

ついでにいうと、OCaml なら破壊的代入ですぐ解決なんですけどね。
まあ関数言語的な考え方はそのぶん習得できましたので、よしとしたいところ。

Haskell の感想

以下、 Haskell の感想

– 状態や逐次処理性を持つアルゴリズム、たとえばシンボルテーブルの管理とかがとっても書きづらい。

State モナドを使っているが、局所的に State モナドを持たせるのができなくて、
関数の上の方まで State モナドの型を”感染”させなくてはならない。
なので結局 do ばっかりのコードになる。

とはいえ、私が単に関数言語的な考え方ができていないので、
Haskell 書きづらいと感じているだけなのかもしれません。
関数型脳な人ならシンボルテーブル管理ももっとじつはうまい方法を思いつくのやも。

C などの手続き言語から入った手続き脳な私は、
どうしてもアルゴリズムを手続きで考えてしまって、
それを関数型にマッピングするという思考になってしまいます。
関数型脳で物ごとを考えるチカラがほしいです。

– コンパイルが遅い

C++ コンパイラ並みにコンパイルが遅い(GHC コンパイラの場合)…

– やっぱコンパイラ書くなら OCaml のほうが書きやすいかも。

だいたいコンパイラは逐次処理が主になるので、遅延評価がかえって邪魔になる。
しかも最近は python や ruby でも関数型的なコードが書けるので、
軽い DSL 向けコンパイラだったら python や ruby を使うのが実は一番手っ取り早いかも。

– cabal(configure みたいなの), darcs(svn みたいなの) のツールチェインで開発は魅力的

とはいえ、Haskell にはいいところもあって、OCaml よりはコードがきれいに書ける(気がする)のと、

http://haskell.org/haskellwiki/How_to_write_a_Haskell_program

にあるように、Cabal でお手軽にプロジェクトのパッケージング&ビルド環境が作れて、
darcs で個人でもお手軽にリビジョン管理できて、
QuickCheck でお手軽にランダムテストが出来て、
haddock でソースコードドキュメンテーションができる。
というよく出来た開発ツールチェインが使えるのはよい。
これだけでも魅力的で Haskell を使おうという理由になる。
# darcs はパッチが多くなると処理がもっさりするのが難点ですが…
# まあその場合は mercurial で置き換えればいいかな.

シェーダ言語も関数型で?

いままでまわりのグラフィックス野郎は C/C++ とかの手続き型言語を利用していて、
関数型言語を知っている人が皆無だったのもあり、関数型言語は私にとって新鮮でした。

何を書きたいかにもよりますが、関数型のほうがよりすっきりとコードが書ける気がします。
(基本的にアルゴリズムそのものが書けるのがよい)

今のシェーディングアルゴリズムやシェーダ言語は、
関数型にすることですっきり書けるような場面が多いのではないかという気がしています。

もちろん、関数型は手続き型にくらべて実行時のパフォーマンスに懸念がある、
という思う方が多いかと思いますが、
これは手続き型から入った人間の先入観でしかないのかもしれません。

DSL(Domain Specific Language)アプローチで適用対象に制限をもうけて、
グラフィックス向けの関数型言語なら、それ向けのコンパイラが頑張ればパフォーマンスには
問題がなくなるのではないかと思っています。

まあグラフィックス用言語に関数型がありえるとしても、
しばらくは python みたいに基本手続き型だけど、部分的に高階関数やλ式などの
関数型な言語仕様を取り込む、みたいなやり方からでしょうかね。

[1] Idea: MUDA, MUltiple Data Accelerator language for high performance computing
http://lucille.atso-net.jp/blog/?p=322

[2] Load to MUDA. SIMD code generation, domain specific language, automatic optimization, functional programming, etc.
http://lucille.atso-net.jp/blog/?p=323

[3] It’s time to move to Haskell to code MUDA?
http://lucille.atso-net.jp/blog/?p=336

Advertisements