Raspberry Pi メモ (17) 各種言語のベンチマーク

Raspberry Pi でプログラミング？

Raspberry Pi は現在(2013/02/02)のところ、XのドライバーがGPUでアクセラレーションされていないため、ウィンドウの消去、移動、スクロールなどでCPUに多くの負荷がかかります。その結果、GUIで使用すると遅く感じます。一方、コンソールに直接OpenGL ES2の出力を表示させると、同じRaspberry Pi とは思えないほどのスピードで3Dアニメーションができます。

3Dのアニメーションを行うゲームなどのアプリケーションは、CやC++でOpenGL ES2のプログラムを作成すれば最も良いパフォーマンスを出せます。しかしプログラミングの学習用として C や C++ を題材にするのは無理があります。たとえば文字列を扱うだけでメモリアドレスやポインタの概念が必要になってきて大変です。 Raspberry Pi のターゲットユーザである子供が皆、CやC++のプログラミングに興味を持ってくれるとも思えません。

Raspberry Pi には、どんな言語が適当なんでしょうか? 公式のSDカードイメージの Raspbian には最初から多くの言語が、すぐ使用できるように入っています。私自身は、最近、Python、C++、LuaJITを使ってOpenGL ES2を試していましたが、いまのところ Raspberry Pi では速度とプログラミングの難易度のバランスで LuaJIT が最も適していると思っています。

今回はまず、OpenGL ES2のプログラムで必要となる 4x4行列の乗算をいろいろな言語で自己中心的にベンチマークをしてみました。あくまで対象を3Dのリアルタイムアニメーションのプログラミング入門としています。今回使った言語は、公式のSDカードイメージの Raspbian(2012-12-16) にすべて最初から含まれています。

結果

まずは結果を示します。LuaJIT は Python の100倍のスピードです。Quaternionと行列の相互変換というような細かい計算を多く行うアプリケーションには非常に有利な処理系です。LuaJIT は ffiモジュールを使うことで、Cのライブラリを自由に呼び出すことができるため、他のアプリケーションに言語機能を提供するための組み込み用言語のLuaとは違って、高速な汎用言語となっています。CPU が遅めの Raspberry Pi にとっては非常に強力なツールになります。

言語	備考	実行時間(秒)
Python 2.7.3	リスト	169.6
Python 2.7.3	NumPy	48.6
Python 3.2.3	リスト	177.4
Python 3.2.3	NumPy	58.9
Lua 5.1.5	-	21.9
LuaJIT 2.0.0	テーブル(JIT)	2.42
	テーブル(JIT OFF)	6.65
	ffi double (JIT)	1.37
	ffi double (JIT OFF)	57.3
gcc 4.6.3	-O0	2.26
gcc 4.6.3	-O2	0.66
アセンブリ	ベクトル演算	0.63

Python

もともとは、教育用に Python を使うために Raspberry Pi の開発が始められたようです。「ラズペリーパイ」の「パイ」は「パイソン」の「パイ」のつもりだったそうです。Python にはctypes という C用の共有ライブラリ内の関数を動的にリンクして呼び出す事のできる強力な機能があります。したがって、Raspberry Pi 固有のlibbcm_host.so などに含まれる関数も ctypes をインポートして、「bcm = ct.CDLL('libbcm_host.so')」のように指定するだけで使えるようになります。

Raspberry Pi には以下のように、Python 2.7.3 と Python 3.2.3 の両方が入っています。

$ python
Python 2.7.3 (default, Jan 13 2013, 11:20:46) 
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>

$ python3
Python 3.2.3 (default, Jul  6 2012, 13:39:51) 
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>

Python のリストを使用

4x4の正方行列を1次元配列で表して計算に使うことにします。Python のデータ型のうち、リストは C の配列とほぼ同じものです。4x4行列の乗算を100万回実行して時間を計っています。

コード

#!/usr/bin/python
# -*- coding: utf-8 -*-

#-------------------------------
# mulmatrix.py
#   2012/02/02 - Jun Mizutani
#-------------------------------
import time

def mulMatrix(mr, ma, mb):
  a0 =ma[ 0];  a1 =ma[ 1];  a2 =ma[ 2];  a3 =ma[ 3]
  a4 =ma[ 4];  a5 =ma[ 5];  a6 =ma[ 6];  a7 =ma[ 7]
  a8 =ma[ 8];  a9 =ma[ 9];  a10=ma[10];  a11=ma[11]
  a12=ma[12];  a13=ma[13];  a14=ma[14];  a15=ma[15]
  b0 =mb[ 0];  b1 =mb[ 1];  b2 =mb[ 2];  b3 =mb[ 3]
  b4 =mb[ 4];  b5 =mb[ 5];  b6 =mb[ 6];  b7 =mb[ 7]
  b8 =mb[ 8];  b9 =mb[ 9];  b10=mb[10];  b11=mb[11]
  b12=mb[12];  b13=mb[13];  b14=mb[14];  b15=mb[15]

  mr[ 0] = a0 * b0 + a4 * b1 +  a8 * b2 + a12 * b3
  mr[ 1] = a1 * b0 + a5 * b1 +  a9 * b2 + a13 * b3
  mr[ 2] = a2 * b0 + a6 * b1 + a10 * b2 + a14 * b3
  mr[ 3] = a3 * b0 + a7 * b1 + a11 * b2 + a15 * b3
  mr[ 4] = a0 * b4 + a4 * b5 +  a8 * b6 + a12 * b7
  mr[ 5] = a1 * b4 + a5 * b5 +  a9 * b6 + a13 * b7
  mr[ 6] = a2 * b4 + a6 * b5 + a10 * b6 + a14 * b7
  mr[ 7] = a3 * b4 + a7 * b5 + a11 * b6 + a15 * b7
  mr[ 8] = a0 * b8 + a4 * b9 +  a8 * b10+ a12 * b11
  mr[ 9] = a1 * b8 + a5 * b9 +  a9 * b10+ a13 * b11
  mr[10] = a2 * b8 + a6 * b9 + a10 * b10+ a14 * b11
  mr[11] = a3 * b8 + a7 * b9 + a11 * b10+ a15 * b11
  mr[12] = a0 * b12+ a4 * b13+  a8 * b14+ a12 * b15
  mr[13] = a1 * b12+ a5 * b13+  a9 * b14+ a13 * b15
  mr[14] = a2 * b12+ a6 * b13+ a10 * b14+ a14 * b15
  mr[15] = a3 * b12+ a7 * b13+ a11 * b14+ a15 * b15

a = [  0.,  1.,  2.,  3., 4.,  5.,  6.,  7.,
       8.,  9., 10., 11., 12., 13., 14., 15.]
b = [  1.,  2.,  3., 4.,  1.,  2.,  3., 4.,
       1.,  2.,  3., 4.,  1.,  2.,  3., 4. ]
c = [  0., 0., 0., 0., 0., 0., 0., 0., 
       0., 0., 0., 0., 0., 0., 0., 0.]

start = time.time()

for i in range(1000000):
    mulMatrix(c, b, a)

print("%f sec\n" % (time.time() - start))
print(c)

実行結果 (Python 2.7.3)

$ time python mulmatrix.py
169.341166 sec

[6.0, 12.0, 18.0, 24.0, 22.0, 44.0, 66.0, 88.0, 38.0, 76.0, 114.0, 152.0, 54.0, 108.0, 162.0, 216.0]

real  2m49.667s
user  2m48.620s
sys 0m0.110s

実行結果 (Python 3.2.3)

$ time python3 mulmatrix.py
176.241442 sec

[6.0, 12.0, 18.0, 24.0, 22.0, 44.0, 66.0, 88.0, 38.0, 76.0, 114.0, 152.0, 54.0, 108.0, 162.0, 216.0]

real  2m57.384s
user  2m57.050s
sys 0m0.110s

他の言語と比べて遅いですが、「100万回の4x4行列の乗算に3分は本当に遅いのか?」という問題です。3次元のアプリケーションでは、物体が運動して、視点が移動する場合には少なくとも2回の乗算が必要になります。「100万 ÷ 180秒 ÷ 秒間60回 ÷ 乗算2回 = 46」なので、行列演算だけで独立した物体が46個が限界になります。色々な最適化や、シェーダの使い方で工夫する余地はありますが、あと10倍ぐらいの速度が欲しい感じです。

NumPy を使用

Python用の高速数値計算モジュールである NumPy を使ってみました。

コード

#!/usr/bin/python
# -*- coding: utf-8 -*-

#-------------------------------
# mulNumPy.py
#   2012/02/02 - Jun Mizutani
#-------------------------------
import time
import numpy as np

a = np.array(
    [[  0.,  1.,  2.,  3.],
     [  4.,  5.,  6.,  7.],
     [  8.,  9., 10., 11.],
     [ 12., 13., 14., 15.]])

b = np.array(
    [[  1.,  2.,  3.,  4.],
     [  1.,  2.,  3.,  4.],
     [  1.,  2.,  3.,  4.],
     [  1.,  2.,  3.,  4.]])

start = time.time()

for i in range(1000000):
    c = np.dot(a, b)

print("%f sec\n" % (time.time() - start))
print(c)

実行結果 (Python 2.7.3)

ench $ time python mulNumPy.py 
46.930748 sec

[[   6.   12.   18.   24.]
 [  22.   44.   66.   88.]
 [  38.   76.  114.  152.]
 [  54.  108.  162.  216.]]

real  0m48.631s
user  0m48.270s
sys 0m0.230s

実行結果 (Python 3.2.3)

$ time python3 mulNumPy.py 
55.821260 sec

[[   6.   12.   18.   24.]
 [  22.   44.   66.   88.]
 [  38.   76.  114.  152.]
 [  54.  108.  162.  216.]]

real  0m58.917s
user  0m58.650s
sys 0m0.170s

Pythonのリストを使った場合より4倍ほど速くなりますが、NumPyのロードに2秒ほどかかっている感じがします。

Lua

Lua は主にアプリケーションやゲームに組み込んで使われる (ゲームコーディング・コンプリート ) ことの多い、高速でコンパクトなスクリプト言語です。最近のプログラミング言語の最大公約数というか、 MECE(ミーシー : mutually exclusive and collectively exhaustive) というか、機能的に必要とされる最小限を取り出したような言語です。Luaは言語仕様が小さく、組み込みの関数や標準ライブラリが多くないため、習得が容易な言語です。

コード

ここでは、Lua のテーブルというデータ型を使って行列を表します。テーブルは連想配列と普通の配列のどちらにも使えます。Lua単独では秒単位の時間計測しかできないので、コマンドの実行時間を計測する bash の組み込みコマンド time を使いました。

#!/usr/bin/lua

-- -----------------------------
-- mulmatrix.lua
--   2012/02/02 - Jun Mizutani
-- -----------------------------

function mulArrayMatrix(mr, ma, mb)
  local a0 =ma[ 1]; local a1 =ma[ 2]; local a2 =ma[ 3]; local a3 =ma[ 4]
  local a4 =ma[ 5]; local a5 =ma[ 6]; local a6 =ma[ 7]; local a7 =ma[ 8]
  local a8 =ma[ 9]; local a9 =ma[10]; local a10=ma[11]; local a11=ma[12]
  local a12=ma[13]; local a13=ma[14]; local a14=ma[15]; local a15=ma[16]
  local b0 =mb[ 1]; local b1 =mb[ 2]; local b2 =mb[ 3]; local b3 =mb[ 4]
  local b4 =mb[ 5]; local b5 =mb[ 6]; local b6 =mb[ 7]; local b7 =mb[ 8]
  local b8 =mb[ 9]; local b9 =mb[10]; local b10=mb[11]; local b11=mb[12]
  local b12=mb[13]; local b13=mb[14]; local b14=mb[15]; local b15=mb[16]

  mr[ 1] = a0 * b0 + a4 * b1 +  a8 * b2 + a12 * b3
  mr[ 2] = a1 * b0 + a5 * b1 +  a9 * b2 + a13 * b3
  mr[ 3] = a2 * b0 + a6 * b1 + a10 * b2 + a14 * b3
  mr[ 4] = a3 * b0 + a7 * b1 + a11 * b2 + a15 * b3
  mr[ 5] = a0 * b4 + a4 * b5 +  a8 * b6 + a12 * b7
  mr[ 6] = a1 * b4 + a5 * b5 +  a9 * b6 + a13 * b7
  mr[ 7] = a2 * b4 + a6 * b5 + a10 * b6 + a14 * b7
  mr[ 8] = a3 * b4 + a7 * b5 + a11 * b6 + a15 * b7
  mr[ 9] = a0 * b8 + a4 * b9 +  a8 * b10+ a12 * b11
  mr[10] = a1 * b8 + a5 * b9 +  a9 * b10+ a13 * b11
  mr[11] = a2 * b8 + a6 * b9 + a10 * b10+ a14 * b11
  mr[12] = a3 * b8 + a7 * b9 + a11 * b10+ a15 * b11
  mr[13] = a0 * b12+ a4 * b13+  a8 * b14+ a12 * b15
  mr[14] = a1 * b12+ a5 * b13+  a9 * b14+ a13 * b15
  mr[15] = a2 * b12+ a6 * b13+ a10 * b14+ a14 * b15
  mr[16] = a3 * b12+ a7 * b13+ a11 * b14+ a15 * b15
end

function printArrayMatrix4(ma)
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[1], ma[5], ma[9],  ma[13]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[2], ma[6], ma[10],  ma[14]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[3], ma[7], ma[11], ma[15]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[4], ma[8], ma[12], ma[16]))
end

local A = { 1, 2, 3, 4,  1, 2, 3, 4,  1, 2, 3, 4,  1, 2, 3, 4}
local B = { 0, 1, 2, 3,  4, 5, 6, 7,  8, 9, 10, 11,  12, 13, 14, 15}
local C = { 0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0}

for i=1,1000000 do
  mulArrayMatrix(C, A, B)
end

printArrayMatrix4(C)

実行

$ time lua mulmatrix.lua
 6.00000e+00  2.20000e+01  3.80000e+01  5.40000e+01
 1.20000e+01  4.40000e+01  7.60000e+01  1.08000e+02
 1.80000e+01  6.60000e+01  1.14000e+02  1.62000e+02
 2.40000e+01  8.80000e+01  1.52000e+02  2.16000e+02

real  0m21.940s
user  0m21.900s
sys 0m0.000s

LuaJIT

LuaJIT は Lua 5.1.5 に対応した JITコンパイラ処理系です。実行時にネイティブコードにコンパイルします。Raspbianのディスクイメージ(2012-12-16-wheezy-raspbian.zip)に入っているのは LuaJIT 2.0.0 直前の LuaJIT 2.0.0-beta11 (2012-10-16) で、初めてRaspberry Piが採用している armhf に対応したバージョンです。インタプリタの手軽さとコンパイラの高速性を備えた非常に高速でコンパクトな処理系です。バイナリのサイズは365キロバイトでPythonの7分の1です。

LuaJIT には標準で FFI (Foreign Function Interface) モジュールという拡張機能を持っていて、その中にCのパーサを持っているため Cのライブラリの宣言(関数や構造体)をほとんどそのままコード中で使うことができます。実行時には共有ライブラリを /etc/ld.so.conf から検索してきて(たぶん)、LuaJITが実行しているコードにリンクできるようです。結果として Python の ctypes と同じことができます。

Luaのテーブルを使う

Luaのテーブルを1次元配列として使います。Luaのテーブルは添字が1から始まるところに注意が必要です。

コード

基本的に上の Luaで書いたコードと同じですが、ここでは内部でも gettimeofday を使って時間を計測しています。これには FFI モジュールを使って C のライブラリの関数を直接使ってPython の time.time() と同じことを実現しています。

#!/usr/bin/luajit

-- -----------------------------
-- mulmatrix1.lua
--   2012/02/02 - Jun Mizutani
-- -----------------------------
ffi = require "ffi"

ffi.cdef[[
    typedef struct timeval {
        long tv_sec;
        long tv_usec;
    } timeval;
    int gettimeofday(struct timeval *tv, void *tz);
]]

function gettimeofday()
  local t = ffi.new("timeval")
  ffi.C.gettimeofday(t, nil)
  return t.tv_sec, t.tv_usec
end

function mulArrayMatrix(mr, ma, mb)
  local a0 =ma[ 1]; local a1 =ma[ 2]; local a2 =ma[ 3]; local a3 =ma[ 4]
  local a4 =ma[ 5]; local a5 =ma[ 6]; local a6 =ma[ 7]; local a7 =ma[ 8]
  local a8 =ma[ 9]; local a9 =ma[10]; local a10=ma[11]; local a11=ma[12]
  local a12=ma[13]; local a13=ma[14]; local a14=ma[15]; local a15=ma[16]
  local b0 =mb[ 1]; local b1 =mb[ 2]; local b2 =mb[ 3]; local b3 =mb[ 4]
  local b4 =mb[ 5]; local b5 =mb[ 6]; local b6 =mb[ 7]; local b7 =mb[ 8]
  local b8 =mb[ 9]; local b9 =mb[10]; local b10=mb[11]; local b11=mb[12]
  local b12=mb[13]; local b13=mb[14]; local b14=mb[15]; local b15=mb[16]

  mr[ 1] = a0 * b0 + a4 * b1 +  a8 * b2 + a12 * b3
  mr[ 2] = a1 * b0 + a5 * b1 +  a9 * b2 + a13 * b3
  mr[ 3] = a2 * b0 + a6 * b1 + a10 * b2 + a14 * b3
  mr[ 4] = a3 * b0 + a7 * b1 + a11 * b2 + a15 * b3
  mr[ 5] = a0 * b4 + a4 * b5 +  a8 * b6 + a12 * b7
  mr[ 6] = a1 * b4 + a5 * b5 +  a9 * b6 + a13 * b7
  mr[ 7] = a2 * b4 + a6 * b5 + a10 * b6 + a14 * b7
  mr[ 8] = a3 * b4 + a7 * b5 + a11 * b6 + a15 * b7
  mr[ 9] = a0 * b8 + a4 * b9 +  a8 * b10+ a12 * b11
  mr[10] = a1 * b8 + a5 * b9 +  a9 * b10+ a13 * b11
  mr[11] = a2 * b8 + a6 * b9 + a10 * b10+ a14 * b11
  mr[12] = a3 * b8 + a7 * b9 + a11 * b10+ a15 * b11
  mr[13] = a0 * b12+ a4 * b13+  a8 * b14+ a12 * b15
  mr[14] = a1 * b12+ a5 * b13+  a9 * b14+ a13 * b15
  mr[15] = a2 * b12+ a6 * b13+ a10 * b14+ a14 * b15
  mr[16] = a3 * b12+ a7 * b13+ a11 * b14+ a15 * b15
end

function printArrayMatrix4(ma)
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[1], ma[5], ma[9],  ma[13]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[2], ma[6], ma[10],  ma[14]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[3], ma[7], ma[11], ma[15]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[4], ma[8], ma[12], ma[16]))
end

local A = { 1, 2, 3, 4,  1, 2, 3, 4,  1, 2, 3, 4,  1, 2, 3, 4}
local B = { 0, 1, 2, 3,  4, 5, 6, 7,  8, 9, 10, 11,  12, 13, 14, 15}
local C = { 0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0 }

local start_sec, start_usec = gettimeofday()

for i=1,1000000 do
  mulArrayMatrix(C, A, B)
end

local end_sec, end_usec = gettimeofday()
local sec = end_sec - start_sec
local usec = end_usec - start_usec
if usec < 0 then 
    usec = usec + 1000000
    sec = sec - 1
end
print(string.format("luajit + array : %f sec", sec+usec/1000000))

printArrayMatrix4(C)

実行

$ time luajit -joff mulmatrix1.lua
luajit + array : 6.622208 sec
 6.00000e+00  2.20000e+01  3.80000e+01  5.40000e+01
 1.20000e+01  4.40000e+01  7.60000e+01  1.08000e+02
 1.80000e+01  6.60000e+01  1.14000e+02  1.62000e+02
 2.40000e+01  8.80000e+01  1.52000e+02  2.16000e+02

real  0m6.646s
user  0m6.610s
sys 0m0.010s

$ time luajit mulmatrix1.lua
luajit + array : 2.394701 sec
 6.00000e+00  2.20000e+01  3.80000e+01  5.40000e+01
 1.20000e+01  4.40000e+01  7.60000e+01  1.08000e+02
 1.80000e+01  6.60000e+01  1.14000e+02  1.62000e+02
 2.40000e+01  8.80000e+01  1.52000e+02  2.16000e+02

real  0m2.418s
user  0m2.390s
sys 0m0.010s

ffi で C の double型の配列を使う

LuaJIT の ffi はCのデータ構造をかなり自由に使うことができます。このバージョンでは行列を倍精度浮動小数点数の1次元配列として表しています。OpenGL ES2 のシェーダに渡す必要のある単精度浮動小数点数の配列も扱えます。

コード

Cの形式の倍精度浮動小数点数の配列なので、インデックス (添字) は 0 から始まります。

#!/usr/bin/luajit

-- -----------------------------
-- mulmatrix2.lua
--   2012/02/02 - Jun Mizutani
-- -----------------------------
ffi = require "ffi"

ffi.cdef[[
    typedef struct timeval {
        long tv_sec;
        long tv_usec;
    } timeval;
    int gettimeofday(struct timeval *tv, void *tz);
]]

function gettimeofday()
  local t = ffi.new("timeval")
  ffi.C.gettimeofday(t, nil)
  return t.tv_sec, t.tv_usec
end

function mulMatrix(mr, ma, mb)
  local a0 =ma[ 0]; local a1 =ma[ 1]; local a2 =ma[ 2]; local a3 =ma[ 3]
  local a4 =ma[ 4]; local a5 =ma[ 5]; local a6 =ma[ 6]; local a7 =ma[ 7]
  local a8 =ma[ 8]; local a9 =ma[ 9]; local a10=ma[10]; local a11=ma[11]
  local a12=ma[12]; local a13=ma[13]; local a14=ma[14]; local a15=ma[15]
  local b0 =mb[ 0]; local b1 =mb[ 1]; local b2 =mb[ 2]; local b3 =mb[ 3]
  local b4 =mb[ 4]; local b5 =mb[ 5]; local b6 =mb[ 6]; local b7 =mb[ 7]
  local b8 =mb[ 8]; local b9 =mb[ 9]; local b10=mb[10]; local b11=mb[11]
  local b12=mb[12]; local b13=mb[13]; local b14=mb[14]; local b15=mb[15]

  mr[ 0] = a0 * b0 + a4 * b1 +  a8 * b2 + a12 * b3
  mr[ 1] = a1 * b0 + a5 * b1 +  a9 * b2 + a13 * b3
  mr[ 2] = a2 * b0 + a6 * b1 + a10 * b2 + a14 * b3
  mr[ 3] = a3 * b0 + a7 * b1 + a11 * b2 + a15 * b3
  mr[ 4] = a0 * b4 + a4 * b5 +  a8 * b6 + a12 * b7
  mr[ 5] = a1 * b4 + a5 * b5 +  a9 * b6 + a13 * b7
  mr[ 6] = a2 * b4 + a6 * b5 + a10 * b6 + a14 * b7
  mr[ 7] = a3 * b4 + a7 * b5 + a11 * b6 + a15 * b7
  mr[ 8] = a0 * b8 + a4 * b9 +  a8 * b10+ a12 * b11
  mr[ 9] = a1 * b8 + a5 * b9 +  a9 * b10+ a13 * b11
  mr[10] = a2 * b8 + a6 * b9 + a10 * b10+ a14 * b11
  mr[11] = a3 * b8 + a7 * b9 + a11 * b10+ a15 * b11
  mr[12] = a0 * b12+ a4 * b13+  a8 * b14+ a12 * b15
  mr[13] = a1 * b12+ a5 * b13+  a9 * b14+ a13 * b15
  mr[14] = a2 * b12+ a6 * b13+ a10 * b14+ a14 * b15
  mr[15] = a3 * b12+ a7 * b13+ a11 * b14+ a15 * b15
end

function printMatrix4(ma)
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[0], ma[4], ma[8],  ma[12]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[1], ma[5], ma[9],  ma[13]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[2], ma[6], ma[10], ma[14]))
  print( string.format(" %10.5e  %10.5e  %10.5e  %10.5e",
           ma[3], ma[7], ma[11], ma[15]))
end

local a = ffi.new("double[?]", 16, 1, 2, 3, 4,  1, 2, 3, 4,  
             1, 2, 3, 4,  1, 2, 3, 4)
local b = ffi.new("double[?]", 16, 0, 1, 2, 3,  4, 5, 6, 7,
                             8, 9, 10, 11,  12, 13, 14, 15)
local c = ffi.new("double[?]", 16, 0, 0, 0, 0,  0, 0, 0, 0,  
                             0, 0, 0, 0,  0, 0, 0, 0)

local start_sec, start_usec = gettimeofday()

for i=1,1000000 do
  mulMatrix(c, a, b)
end

local end_sec, end_usec = gettimeofday()
local sec = end_sec - start_sec
local usec = end_usec - start_usec
if usec < 0 then 
    usec = usec + 1000000
    sec = sec - 1
end
print(string.format("luajit + ffi : %f sec", sec+usec/1000000))

printMatrix4(c)

jit off 実行

$ time luajit -joff mulmatrix2.lua
luajit + ffi : 57.263448 sec
 6.00000e+00  2.20000e+01  3.80000e+01  5.40000e+01
 1.20000e+01  4.40000e+01  7.60000e+01  1.08000e+02
 1.80000e+01  6.60000e+01  1.14000e+02  1.62000e+02
 2.40000e+01  8.80000e+01  1.52000e+02  2.16000e+02

real  0m57.287s
user  0m57.180s
sys 0m0.030s

jit on 実行

$ time luajit mulmatrix2.lua
luajit + ffi : 1.343749 sec
 6.00000e+00  2.20000e+01  3.80000e+01  5.40000e+01
 1.20000e+01  4.40000e+01  7.60000e+01  1.08000e+02
 1.80000e+01  6.60000e+01  1.14000e+02  1.62000e+02
 2.40000e+01  8.80000e+01  1.52000e+02  2.16000e+02

real  0m1.368s
user  0m1.340s
sys 0m0.010s

C とアセンブリ

Raspberry Pi メモ (05) とほぼ同じコードですが、比較のためにここでもリストを載せておきます。

C の部分

アセンブリで書いたオブジェクトファイルをリンクしている以外は、他の言語と同じ事をしています。

コード

行列が配列ではなく、構造体にしています。

//--------------------------------------------------------
// Multiply 4x4 Matrix
//   2013/02/02 Jun Mizutani
// as -mcpu=arm1176jzf-s -mfpu=vfp mat4x4.s -o mat4x4asm.o
// gcc -Wall -o mat4 mat4.c mat4x4asm.o
// gcc -O2 -Wall -o mat4-O2 mat4.c mat4x4asm.o
//--------------------------------------------------------

#include <stdio.h>
#include <sys/time.h>

typedef struct Matrix4 {
    double m0,m1,m2,m3, m4,m5,m6,m7, m8,m9,m10,m11, m12,m13,m14,m15;
} Matrix4;

Matrix4 a = { 1, 2, 3, 4,  1, 2, 3, 4,  1, 2, 3, 4,     1, 2, 3, 4 };
Matrix4 b = { 0, 1, 2, 3,  4, 5, 6, 7,  8, 9, 10, 11,  12, 13, 14, 15 };
Matrix4 c = { 0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,     0, 0, 0, 0 };

extern void MulMatrix4x4asm(Matrix4 *result, Matrix4 *a, Matrix4 *b);

void MulMatrix4(Matrix4 *a, Matrix4 *b, Matrix4 *c)
{
    c->m0 = a->m0 * b->m0 + a->m4 * b->m1 + a->m8 * b->m2 + a->m12 * b->m3;
    c->m1 = a->m1 * b->m0 + a->m5 * b->m1 + a->m9 * b->m2 + a->m13 * b->m3;
    c->m2 = a->m2 * b->m0 + a->m6 * b->m1 + a->m10* b->m2 + a->m14 * b->m3;
    c->m3 = a->m3 * b->m0 + a->m7 * b->m1 + a->m11* b->m2 + a->m15 * b->m3;

    c->m4 = a->m0 * b->m4 + a->m4 * b->m5 + a->m8 * b->m6 + a->m12 * b->m7;
    c->m5 = a->m1 * b->m4 + a->m5 * b->m5 + a->m9 * b->m6 + a->m13 * b->m7;
    c->m6 = a->m2 * b->m4 + a->m6 * b->m5 + a->m10* b->m6 + a->m14 * b->m7;
    c->m7 = a->m3 * b->m4 + a->m7 * b->m5 + a->m11* b->m6 + a->m15 * b->m7;

    c->m8 = a->m0 * b->m8 + a->m4 * b->m9 + a->m8 * b->m10+ a->m12 * b->m11;
    c->m9 = a->m1 * b->m8 + a->m5 * b->m9 + a->m9 * b->m10+ a->m13 * b->m11;
    c->m10= a->m2 * b->m8 + a->m6 * b->m9 + a->m10* b->m10+ a->m14 * b->m11;
    c->m11= a->m3 * b->m8 + a->m7 * b->m9 + a->m11* b->m10+ a->m15 * b->m11;

    c->m12= a->m0 * b->m12+ a->m4 * b->m13+ a->m8 * b->m14+ a->m12 * b->m15;
    c->m13= a->m1 * b->m12+ a->m5 * b->m13+ a->m9 * b->m14+ a->m13 * b->m15;
    c->m14= a->m2 * b->m12+ a->m6 * b->m13+ a->m10* b->m14+ a->m14 * b->m15;
    c->m15= a->m3 * b->m12+ a->m7 * b->m13+ a->m11* b->m14+ a->m15 * b->m15;
}

void PrintMatrix4(Matrix4 *a)
{
    printf("%10.5e %10.5e %10.5e %10.5e \n", a->m0, a->m4, a->m8,  a->m12);
    printf("%10.5e %10.5e %10.5e %10.5e \n", a->m1, a->m5, a->m9,  a->m13);
    printf("%10.5e %10.5e %10.5e %10.5e \n", a->m2, a->m6, a->m10, a->m14);
    printf("%10.5e %10.5e %10.5e %10.5e \n\n", a->m3, a->m7, a->m11, a->m15);
}

int main(int argc, char* argv[])
{
  struct timeval start, now;
  struct timezone tzone;
  double elapsedTime = 0.0;
  int i;

  printf("start C version\n");
  gettimeofday(&start, &tzone);

  for(i=0; i<1000000; i++) MulMatrix4(&a, &b, &c);

  gettimeofday(&now, &tzone);
  elapsedTime = (double)(now.tv_sec - start.tv_sec) +
    (double)(now.tv_usec - start.tv_usec)/1000000.0;
  printf("%4.6f sec\n", elapsedTime);

  PrintMatrix4(&c);

  printf("start ASM version\n");
  gettimeofday(&start, &tzone);

  for(i=0; i<10000000; i++) MulMatrix4x4asm(&c, &a, &b);

  gettimeofday(&now, &tzone);
  elapsedTime = (double)(now.tv_sec - start.tv_sec) +
    (double)(now.tv_usec - start.tv_usec)/1000000.0;
  printf("%4.6f sec\n", elapsedTime);

  PrintMatrix4(&c);

  return 0;
}

アセンブリ部分

VFPv2のベクトルモードを使って行列の乗算を行なっています。詳細は Raspberry Pi メモ (05) を参照して下さい。

コード

@----------------------------------------------------------
@ mat4x4.s   2012/08/19 Jun Mizutani 
@ as -mcpu=arm1176jzf-s -mfpu=vfp mat4x4.s -o mat4x4asm.o
@ gcc -O2 -Wall -o mat4 mat4.c mat4x4asm.o
@---------------------------------------------------------

.text
.global MulMatrix4x4asm

@ MulMatrix4x4asm(Matrix4 *result, Matrix4 *a, Matrix4 *b)
@ [result] = [a] x [b]
@ [  r0  ] = [r1] x [r2]
        .align  3
MulMatrix4x4asm:
        stmfd   sp!, {r0-r3, lr}
        vpush   {d0-d15}
        bl      vector4_mode            @ LEN=4: vector mode
        mov     r3, r0
        fldmiad r2!, {d0-d3}            @ vec_e
        mov     r0, r1                  @ vec_a
        fldmiad r0!, {d8-d11}           @ load a[0_3]
        fmuld   d4, d8, d0              @ w[i] =        a[i] * e0
        fldmiad r0!, {d12-d15}          @ load b[0-3]
        fmacd   d4, d12, d1             @ w[i] = w[i] + b[i] * e1
        fldmiad r0!, {d8-d11}           @ load c[0-3]
        fmacd   d4, d8, d2              @ w[i] = w[i] + c[i] * e2
        fldmiad r0, {d12-d15}           @ load d[0-3]
        fmacd   d4, d12, d3             @ w[i] = w[i] + d[i] * e3
        fstmiad r3!, {d4-d7}            @ store w[i]
        fldmiad r2!, {d0-d3}            @ vec_f
        mov     r0, r1                  @ vec_a
        fldmiad r0!, {d8-d11}
        fmuld   d4, d8, d0              @ x[i] =        a[i] * f0
        fldmiad r0!, {d12-d15}
        fmacd   d4, d12, d1             @ x[i] = x[i] + b[i] * f1
        fldmiad r0!, {d8-d11}
        fmacd   d4, d8, d2              @ x[i] = x[i] + c[i] * f2
        fldmiad r0, {d12-d15}
        fmacd   d4, d12, d3             @ x[i] = x[i] + d[i] * f3
        fstmiad r3!, {d4-d7}            @ store x[i]
        fldmiad r2!, {d0-d3}            @ vec_g
        mov     r0, r1                  @ vec_a
        fldmiad r0!, {d8-d11}
        fmuld   d4, d8, d0              @ y[i] =        a[i] * g0
        fldmiad r0!, {d12-d15}
        fmacd   d4, d12, d1             @ y[i] = y[i] + b[i] * g1
        fldmiad r0!, {d8-d11}
        fmacd   d4, d8, d2              @ y[i] = y[i] + c[i] * g2
        fldmiad r0, {d12-d15}
        fmacd   d4, d12, d3             @ y[i] = y[i] + d[i] * g3
        fstmiad r3!, {d4-d7}            @ store y[i]
        fldmiad r2, {d0-d3}             @ vec_h
        mov     r0, r1                  @ vec_a
        fldmiad r0!, {d8-d11}
        fmuld   d4, d8, d0              @ z[i] =        a[i] * h0
        fldmiad r0!, {d12-d15}
        fmacd   d4, d12, d1             @ z[i] = z[i] + b[i] * h1
        fldmiad r0!, {d8-d11}
        fmacd   d4, d8, d2              @ z[i] = z[i] + c[i] * h2
        fldmiad r0, {d12-d15}
        fmacd   d4, d12, d3             @ z[i] = z[i] + d[i] * h3
        fstmiad r3, {d4-d7}

        bl      scalar_mode             @ LEN=1: scalar mode
        vpop    {d0-d15}
        ldmfd   sp!, {r0-r3, pc}

@-------------------------------------------------------------------------
@ set FPSTR STRIDE[21:20] and LEN[18:16]
@-------------------------------------------------------------------------
scalar_mode:
        stmfd   sp!, {r2, lr}
        fmrx    r2, FPSCR              @ FPSCR --> r10
        bic     r2, r2, #0x00370000    @ clear (STRIDE=1 and LEN=1)
        fmxr    FPSCR, r2              @ r10 --> FPSCR
        ldmfd   sp!, {r2, pc}

ビルド

#!/bin/sh

as -mcpu=arm1176jzf-s -mfpu=vfp mat4x4.s -o mat4x4asm.o
gcc -O0 -Wall -o mulmat mat4.c mat4x4asm.o
gcc -O2 -Wall -o mulmat-O2 mat4.c mat4x4asm.o

実行

jun@raspi512 ~/MatrixBench $ time ./build

real  0m2.792s
user  0m2.360s
sys 0m0.240s
jun@raspi512 ~/MatrixBench $ ./mulmat
start C version
2.261465 sec
6.00000e+00 2.20000e+01 3.80000e+01 5.40000e+01 
1.20000e+01 4.40000e+01 7.60000e+01 1.08000e+02 
1.80000e+01 6.60000e+01 1.14000e+02 1.62000e+02 
2.40000e+01 8.80000e+01 1.52000e+02 2.16000e+02 

start ASM version
0.677662 sec
6.00000e+00 2.20000e+01 3.80000e+01 5.40000e+01 
1.20000e+01 4.40000e+01 7.60000e+01 1.08000e+02 
1.80000e+01 6.60000e+01 1.14000e+02 1.62000e+02 
2.40000e+01 8.80000e+01 1.52000e+02 2.16000e+02 

jun@raspi512 ~/MatrixBench $ ./mulmat-O2 
start C version
0.662011 sec
6.00000e+00 2.20000e+01 3.80000e+01 5.40000e+01 
1.20000e+01 4.40000e+01 7.60000e+01 1.08000e+02 
1.80000e+01 6.60000e+01 1.14000e+02 1.62000e+02 
2.40000e+01 8.80000e+01 1.52000e+02 2.16000e+02 

start ASM version
0.628814 sec
6.00000e+00 2.20000e+01 3.80000e+01 5.40000e+01 
1.20000e+01 4.40000e+01 7.60000e+01 1.08000e+02 
1.80000e+01 6.60000e+01 1.14000e+02 1.62000e+02 
2.40000e+01 8.80000e+01 1.52000e+02 2.16000e+02

Cもアセンブラも同じように最速です。Pythonと比べると300倍近い差になります。かなり頑張って書いたアセンブリのコードとLuaJITが2倍程度しか速くないのはショックです。逆に LuaJIT の速度を信じて安心して使ってもいいとも言えます。

ダウンロード

ここに掲載したコード(MatrixBench.tar.gz 9KB)を置いておきます。

$ tar ztvvf MatrixBench.tar.gz
drwxr-xr-x jun/jun           0 2013-02-02 16:58 MatrixBench/
-rw-r--r-- jun/jun         917 2013-02-02 16:57 MatrixBench/mat4x4asm.o
-rwxr-xr-x jun/jun         156 2013-02-02 14:48 MatrixBench/build
-rw-r--r-- jun/jun        3674 2013-02-02 06:06 MatrixBench/mat4x4.s
-rwxr-xr-x jun/jun        2975 2013-02-02 12:36 MatrixBench/mulmatrix2.lua
-rwxr--r-- jun/jun        1817 2013-02-02 15:02 MatrixBench/mulmatrix.py
-rwxr-xr-x jun/jun        7357 2013-02-02 16:57 MatrixBench/mulmat-O2
-rwxr--r-- jun/jun        3149 2013-02-02 13:09 MatrixBench/mat4.c
-rwxr-xr-x jun/jun        2269 2013-02-02 12:23 MatrixBench/mulmatrix.lua
-rwxr--r-- jun/jun         561 2013-02-02 14:43 MatrixBench/mulNumPy.py
-rwxr-xr-x jun/jun        8749 2013-02-02 16:57 MatrixBench/mulmat
-rwxr-xr-x jun/jun        2853 2013-02-02 12:35 MatrixBench/mulmatrix1.lua

Raspberry Pi(ラズベリーパイ) メモ (17) 各種言語のベンチマーク

Raspberry Pi でプログラミング？

結果

Python

Python のリストを使用

コード

実行結果 (Python 2.7.3)

実行結果 (Python 3.2.3)

NumPy を使用

コード

実行結果 (Python 2.7.3)

実行結果 (Python 3.2.3)

Lua

コード

実行

LuaJIT

Luaのテーブルを使う

コード

実行

ffi で C の double型の配列を使う

コード

jit off 実行

jit on 実行

C とアセンブリ

C の部分

コード

アセンブリ部分

コード

ビルド

実行

ダウンロード

「Minecraft - Raspberry Pi版」に続く...

Jun's Homepage

ラズベリーパイメモ

ARM アセンブリ

このページの目次

Raspberry Pi(ラズベリーパイ) メモ (17) 各種言語のベンチマーク

Raspberry Pi でプログラミング？

結果

Python

Python のリストを使用

コード

実行結果 (Python 2.7.3)

実行結果 (Python 3.2.3)

NumPy を使用

コード

実行結果 (Python 2.7.3)

実行結果 (Python 3.2.3)

Lua

コード

実行

LuaJIT

Luaのテーブルを使う

コード

実行

ffi で C の double型の配列を使う

コード

jit off 実行

jit on 実行

C とアセンブリ

C の部分

コード

アセンブリ部分

コード

ビルド

実行

ダウンロード

「Minecraft - Raspberry Pi版」 に続く...

Jun's Homepage

ラズベリーパイ メモ

ARM アセンブリ

このページの目次

「Minecraft - Raspberry Pi版」に続く...

ラズベリーパイメモ