Syoyo Fujita's Blog

raytracing monte carlo

Category: compiler

no divdi3 in x64 libgcc?

Win64 環境(mingw x64)で LLVM 2.6 をビルドしようとして、
libgcc の奇妙な?構成に気がつきました.

libgcc には、除算器がないハードでも、除算をソフトウェア実現するコードなどがあります.
そこには、たとえば __divdi3 というものがあるのですが、
これがなぜか 64bit libgcc では, object 名は divdi3 なのに、実装の関数は divti3 と、名前が変わっています.

nm で libgcc を調べると、以下のようになります.


_divdi3.o:
...
0000000000000000 T ___divti3

はてなんででしょう?
ちなみに divdi3 は 64bit integer の除算、divti3 は 128bit integer の除算という役割になっています(なっているはず).
なので、divti3 は、object 名も divti3.o にして、divdi3.o には divdi3 の関数を実装するべきだとおもうのですが…

gcc ML にでも聞いてみますかね.

World’s first SIMD ray-triangle intersection with ARM/NEON

I’ve wrote possibly world’s first SIMD ray-triangle intersection code with ARM/NEON instruction.

The code runs finely on iPod touch 3G(ARM Cortex) and iPhone 3GS.
The performance on iPod touch 3G is around 348 cycles per isect4().
# of assembly instruction of isect4() is around 140, thus the CPI = 348/140 = 2.48, which is not an efficient number for me(CPI should be near to 1.0).

Note that there’s still a room to optimize the code.
e.g. use fmad instruction instead of fadd/fmul combination.

Related articles

– SIMD ray-triangle intersection in Intel/AVX
http://lucille.atso-net.jp/blog/?p=649

– SIMD ray-triangle intersection in OpenCL
http://lucille.atso-net.jp/blog/?p=910

[Google groups] SIMD-cafe for SIMD lovers

なぜ SIMD が重要か

プロセッサのクロック数は頭打ち…
マルチコア化ではコアが増えるたびにコア間の通信がネックになる問題がある.
専用プロセッサを作るにも, 固定したアルゴリズム&量が出ないとビジネスにならないので,
やはり CPU 的にソフトなプロセッサが求められる.
一方で、エコ時代であるので省電力が求められる.
そんな HW の制約の中で、SW(アプリ)は富豪的プログラミングで HW の向上に頼るだけでは対応しきれなくなってきています.
特に認識処理、グラフィックス、データベース処理でパフォーマンスを出していくことが求められている.

このような、現在のプロセッサ(半導体)の技術的壁やプロセッサデザインの収益モデル(ビジネスモデル)を考えると、
特に SIMD 演算器とその活用が重要になってきます.

実際, x86/AVX(wider SIMD), iPhone 3GS(ARM/NEON), Atom/SSE が登場してきています.
アプリを SIMD 最適化する必要性は今後また重要になってくるでしょう.
(SIMD 最適化の最初のブームは、 MMX/SSE/SSE2 の出てきた 2000 年頃でしょうか)

しかし、どんなにコンパイラが賢くても、効率的な SIMD コードを自動で吐いてくれるのはまだまだ難しい.

そこで、プログラマがいかにアルゴリズムを SIMD フレンドリーにするか、
SIMD 演算器をいかにまく使ってプログラムを最適化していくかが、
今また大きく求められている、そんな気がしています.

また、ときどき SIMD 野郎と話すことがあるのですが、
結構みんな「ほかに SIMD のことを濃く話すことができるひとが居なくて…」
と言われることもあり, SIMD 野郎のための google groups を作ってみました.

SIMD cafe
http://groups.google.co.jp/group/simd-cafe

SIMD 野郎の参加をお待ちしております.

OpenCL ray-triangle intersection kernel.

openclraycastsc.png

I’ve wrote possibly World’s first ray-triangle intersection kernel in OpenCL.
I am now working with porting full raytracing kernel(traversal, shader, isector, etc) into OpenCL.
This article is just exposing the small part of my professional(?) OpenCL work to you.

Here’s the sample code of ray-triangle intersection kernel in OpenCL.


//
// Copyright 2009 Syoyo Fujita.
//
typedef struct _ray4_t
{
    float4 rox, roy, roz;
    float4 rdx, rdy, rdz;
    float4 t;
} ray4_t;

typedef struct _triangle4_t
{
    float4 v0x, v0y, v0z;
    float4 e1x, e1y, e1z;
    float4 e2x, e2y, e2z;
} triangle4_t;

typedef struct _triangle_t
{
    float v[3][4];  // (x,y,z,w) * 3
} triangle_t;

static inline float4
mycross4(float4 a, float4 b, float4 c, float4 d)
{
    return ((a * c) - (b * d));
}

static inline float4
mydot4(float4 ax, float4 ay, float4 az, float4 bx, float4 by, float4 bz)
{
    return (ax * bx + ay * by + az * bz);
}

static int4
isect4(
    float4 *t_out,
    float4 *u_out,
    float4 *v_out,
    ray4_t ray,
    triangle4_t tri)
{
    const float4 px = mycross4(tri.e2z, tri.e2y, ray.rdy, ray.rdz);
    const float4 py = mycross4(tri.e2x, tri.e2z, ray.rdz, ray.rdx);
    const float4 pz = mycross4(tri.e2y, tri.e2x, ray.rdx, ray.rdy);

    const float4 sx = ray.rox - tri.v0x;
    const float4 sy = ray.roy - tri.v0y;
    const float4 sz = ray.roz - tri.v0z;

    const float4 vone  = (float4)(1.0f);
    const float4 vzero = (float4)(0.0f);
    const float4 veps  = (float4)(1.0e-6f);

    const float4 det = mydot4(px, py, pz, tri.e1x, tri.e1y, tri.e1z);
    const float4 invdet = (float4)(1.0f) / det;


    const float4 qx = mycross4(tri.e1z, tri.e1y, sy, sz);
    const float4 qy = mycross4(tri.e1x, tri.e1z, sz, sx);
    const float4 qz = mycross4(tri.e1y, tri.e1y, sx, sy);

    const float4 u = mydot4(sx, sy, sz, px, py, pz) * invdet;
    const float4 v = mydot4(ray.rdx, ray.rdy, ray.rdz, qx, qy, qz) * invdet;
    const float4 t = mydot4(tri.e2x, tri.e2y, tri.e2z, qx, qy, qz) * invdet;

    int4 mask = (fabs(det) > veps) & ((u+v)  vzero) & (v > vzero)
              & (t > vzero) & (t < ray.t);


    (*t_out) = t;
    (*u_out) = u;
    (*v_out) = v;

    return mask;

}

__kernel
void
main(
    __global uint *out)
{
    int x = get_global_id(0);
    int y = get_global_id(1);
    int w = get_global_size(0);
    int h = get_global_size(1);

    ray4_t ray4;

    ray4.rox = 2.0f * ((x - w) / (float)w) + 1.0f;
    ray4.roy = 2.0f * ((y - h) / (float)h) + 1.0f;
    ray4.roz = (float4)(0.0f);

    ray4.rdx = (float4)(0.0f);
    ray4.rdy = (float4)(0.0f);
    ray4.rdz = (float4)(-1.0f);

    ray4.t = (float4)(1.0e+30f);



    triangle_t tri;
    triangle4_t tri4;

    tri.v[0][0] = -0.5f;
    tri.v[0][1] = -0.5f;
    tri.v[0][2] = -1.0f;
    tri.v[1][0] =  0.1f;
    tri.v[1][1] =  0.75f;
    tri.v[1][2] = -1.0f;
    tri.v[2][0] =  0.5f;
    tri.v[2][1] = -0.75f;
    tri.v[2][2] = -1.0f;

    tri4.v0x = (float4)(tri.v[0][0]);
    tri4.v0y = (float4)(tri.v[0][1]);
    tri4.v0z = (float4)(tri.v[0][2]);
    tri4.e1x = (float4)(tri.v[1][0] - tri.v[0][0]);
    tri4.e1y = (float4)(tri.v[1][1] - tri.v[0][1]);
    tri4.e1z = (float4)(tri.v[1][2] - tri.v[0][2]);
    tri4.e2x = (float4)(tri.v[2][0] - tri.v[0][0]);
    tri4.e2y = (float4)(tri.v[2][1] - tri.v[0][1]);
    tri4.e2z = (float4)(tri.v[2][2] - tri.v[0][2]);

    float4 t, u, v;
    int4 mask = isect4(&t, &u, &v, ray4, tri4);

    //int r = mask.x & 0xff;
    //int g = mask.y & 0xff;

    int r = (int)(t.x * 255);
    int g = (int)(u.x * 255);
    int b = (int)(v.x * 255);

    r = mask.x & r;
    g = mask.x & g;
    b = mask.x & b;

    out[x+y*get_global_size(0)] = (uint)(r | (g<< 8) | (b << 16) | (255<<24));
}

Similar articles

World’s first SIMD ray-triangle intersection in AVX?
http://lucille.atso-net.jp/blog/?p=649

KLEE : Unassisted and Automatic Generation of High-Coverage Tests for Complex Systems Programs

KLEE は LLVM を使ったシンボリック実行ツール.

http://klee.llvm.org/

シンボリック実行とは?

シンボリック実行(Symbolic execution)とは、論文によれば以下のようなシステムだそう.

Instead of running code on manually- or randomly-constructed input,
they run it on symbolic input initially allowed to be “anything.”
They substitute program inputs with symbolic values and replace
corresponding concrete program operations with ones that manipulate symbolic values.

このシステムの使い道として論文が述べているのは主に二つ.
カバレッジツールと、自動バグ発見ツールである.

カバレッジツール

シンボリック実行を行うことにより、KLEE が自動生成するテストデータでは 90% 以上と高いコードカバレッジを実現する.
たとえば Coreutils に対するカバレッジテストでは、KLEE はわずか 89 時間をかけることで、
15 年以上手動で積み上げられてきたテストスイートを上回るカバレッジを達成している.

自動バグ発見ツール

KLEE はバグ発見ツールとしても利用でき、Coreutils に 15 年潜んできた致命的なバグを発見できたという.

まとめ

KLEE およびシンボリック実行アルゴリズム、なかなかよさそうですね.
レンダラの検証にも使えないかなと思います.
たとえばシーンファイルのパーサに適用とか.

ただ、ちゃんとレンダラアプリの内部までのカバレッジをあげるには、coreutils のようなある意味テキスト操作系の
アプリとはまた違った戦略が必要になりますので、KLEE そのままの適用では難しいでしょうね.
そこらへん、新しくロジックを提案できないかなぁと思います.

Ocelot : Binary translator for PTX.

(via gpupapers)

ocelotfigure.png

http://code.google.com/p/gpuocelot/

http://code.google.com/p/gpuocelot/wiki/References

(papers)

CUDA で使われている PTX(中間言語表現)の、バイナリトランスレータ.
名前は Ocelot.

PTX の時点で、プログラムは並列性を生かしたものになっているので、
要は PTX コードを変換するツールを書けば,
高級言語(CUDA)の段階からうまーく並列化されたコードが書けていろんなプラットフォームで動かせるぜ、という発想らしい。

PTX へのバイナリトランスレータを作ることで、

– PTX のスレッド階層はメニーコアにもマップ可能
– 変換のテクニックは、メモリレイテンシの隠蔽にも使うことができる。
– GPU データ構造を効率的にエミュレートしたり似たようなプロセッサにマップできる

なことを示している。
例として、PTX -> SPE(Cell) アセンブラへのバイナリトランスレータを作り、
うまくいったということを示している.

Binary Translator であることの利点

Ocelot は PTX のバイナリトランスレータでり、既存の MCUDA などの CUDA を CPU で動かすなどの手法と比べて以下の 2 点の利点があると述べている.

– スレッドの階層レベルをランタイム時に動的に割り当てることができる.
– 実行時情報を利用して、hot path 最適化などが行える.

気になる点

SPU に変換したときのパフォーマンスが書かれていないですね。
sin とか同期プリミティブなどの PTX(GPU) 専用の命令が SPU には無いので,
これを SPU ではエミュレートする必要があり命令長がそれなりに増えそう.
これがどれだけパフォーマンスに影響を与えるか気になります.

また、浮動小数点の演算精度が PTX(NV) と SPU では微妙に異なるので、
完全一致な結果を得られることができません。
ゲーム用などならいいかもしれませんが、
精度を求められる用途だとどうかなーという気がします。

gpuocelet

Ocelot は open source で公開されており、
LLVM IR をサポートするなど論文では取り上げてないものも含まれている.
将来的には、PTX から x86/LLVM/LRB/OpenCL などへのバイナリ変換ツールを目指しているようだ。

まとめ

中間言語に PTX を使っちゃう手法、なかなかいいんじゃないかと思います。
ユーザとしては CUDA プログラムをするか、CUDA コードを吐く DSL コンパイラを作ることで、いろいろなプロセッサで動く並列コードが書けるようになるわけですし。
CPU に出すことで、デバッグも楽になるんじゃないかな.

また、PTX -> OpenCL という方法も取れるわけで、既存のアプリの OpenCL 化にも使えそうです.
このとき、単にトランスレートするのではなく、動的実行してプロファイリングを取り、ターゲットプロセッサに最適な OpenCL コードを吐いたりとかできそう.
(OpenCL は中間言語の規定がないので、このような動的最適化をするのが難しい)

AMD OpenCL(x86) v.s. MUDA

Recently, AMD released OpenCL drivers which runs on x86 CPU.
I think it’d be better idea to measure performance of Black-Scholes compuation on AMD’s OpenCL(x86) and MUDA(SSE x86).

Both was run on 1.86 GHz Core2 Linux 32bit.

AMD’s OpenCL(x86)


Option samples           Time taken(sec)          Options / sec
10240000                 4.497                    2.27707e+06

AMD’s OpenCL results in 2.3 M ops/sec.

MUDA version

MUDA version of BlackScholes are taken from MUDA’s source code tree.
https://lucille.svn.sourceforge.net/svnroot/lucille/angelina/haskellmuda/examples/BlackScholes/


[Setup] Generating input data ...
[Setup] Generating input data ... DONE
[Measure] Computing reference ...
[Measure] Computing reference ... DONE
[Measure] Computing with MUDA ...
[Measure] Computing with MUDA ... DONE
[Perf] CPU  = 4781.213000 (msec)  4.283432 MOps/sec
[Perf] MUDA = 781.397300 (msec)  26.209458 MOps/sec
L1 norm: 1.853851E-07
Max absolute error: 1.335144E-05
TEST PASSED

MUDA computes B/S at 26 M ops/sec, which is 10x faster than OpenCL(x86) version!

Because MUDA is open source, you could find why MUDA is faster than AMD’s OpenCL code, hehe…

[Ja]

AMD から OpenCL ベータ実装がリリースされ、x86 でもそのカーネルが動くということで、
MUDA(SSE x86 出力) と比較してみました.
比較対象は BlackScholes 計算式.

結果は…

AMD’s OpenCL(x86) 2.3 M ops/sec
MUDA 26 M ops/sec

… なんだ、圧倒的じゃないか、我が MUDA は 😉
んー、なんで AMD の OpenCL はこんなに遅いんですかねー.
サンプル数は 10240000 個と、データ転送や API コールがネックにならないようにしてみたのですが.

ちなみに、AMD の OpenCL SDK の配布物をいろいろ見てみたのですが、
OpenCL 言語のコンパイラに LLVM 使っていますね.

OSL(Open Shading Language)

(via Blender to renderman group)

http://opensource.imageworks.com/

Open Shading Language (OSL) is a small but rich language for programmable shading in advanced renderers and other applications. OSL is similar to C, as well as other shading languages, however, it is specifically designed for advanced rendering algorithms with features such as radiance closures, BRDFs, and deferred ray tracing as first-class concepts.

Imageworks からいくつかの OpenSource プロジェクトが公開され、
そのうちのひとつが OSL(Open Shading Language).

オープンな仕様で、レイトレベースで、コンパイラもオープンソースで提供されるらしい.

オープンなレイトレシェーダ言語はこちらも考えていたので、さて OSL がどれくらい使われるようになるかですね.
よい仕様でこれが普及してくれば、OSL にスイッチして lucille に取り入れるかもしれません.
まー、でも C ベースの言語仕様で、またコンパイラが Haskell で書かれていたりするわけでもないので、
OSL はフツーな言語でフツーな実装になりそうなのがちょっと残念かなぁ.
(ま、実際に仕様やコードが公開されてからいろいろ判断しますか)

MILEPOST GCC: optimizing a compiler by machine learning

http://www.milepost.eu/

機械学習を用いて、いくつかサンプルプログラムを走らせ学習し、そのプロセッサに最適なコンパイルオプションの組み合わせを導出するというもの.
単純に “-O3” 最適化オプションを指定するものと比べて、MILEPOST GCC を使うと 11% ~ 16% 程度の高速化を得ることができる. しかもコンパイラオプションのどの組み合わせが有効かはこのフレームワークが勝手に導出してくれるのでとってもお手軽.

たとえば遺伝的アルゴリズムを利用した ACOVEA のように、

http://www.coyotegulch.com/products/acovea/

このような自動で最適化オプションを探すという手法はありましたが、
探索区間が大きかったり試行に時間がかかるなどでチューニング時間がかかり広く使える手法にはほど遠かったとのこと.
(ACOVEA は確か 1 日近くかかる)

MILEPOST GCC は、その問題を機械学習を使うことで効率化し、
またフレームワークを GCC plugin で実現することで実務的に使えるようにしています.

実際ソースを落としていじろうとしましたが、
最後の最後らへんで、うまく動いてくれませんでした.

ちなみに、学習データは ctuning.org でデータベース化していろいろ面白いことをやっていこうとするようです.

http://ctuning.org/wiki/index.php/Main_Page

Wirbel, statically typed Python-like language.

http://mathias-kettner.de/wirbel_index.html

Python 的な文法だけど、静的型付けな言語.
もちろん(?) モダンな言語ということで型付けは型推論により行われます.

出来るコードはネイティブバイナリなので高速動作.
また内部データ型は C と合わせてあるので, C/C++ との統合も用意.

な、なんかすごくね?…

こんな感じで書いて…


def main():

        a = 1.0
        b = 3.0

        c = [1, 2, 3]

        for i in c:
                print(i)

        print("a + b = %f" % (a+b))

        print("The World! Time stops...")

main()

こんな感じで実行します.


$ wirbel main.w
1
2
3
a + b = 4.000000
The World! Time stops...

$ wic main.w
$ file main
main: Mach-O executable i386
$ ./main
1
2
3
a + b = 4.000000
The World! Time stops...

print がかっこ付きなのは Python3.0 っぽくていいんじゃないでしょうか.

これはなかなか期待できそうな言語ですね.

Wirbel, 次世代シーンファイルフォーマットの記述言語の候補としても考えておきたい.
ライセンスが GPL ってのが気がかりですが.