Raspberry Pi で Java8 プログラミング (2014/12/29)

Raspberry Pi の OS、Raspbian の 2014年9月版以降には、Java の最新版の JDK 8 が最初からディスクイメージに含まれています。 Raspberry Pi だけで最新の Java8 のプログラミングができます。 これまで私は Java には近づかなかったのですが、いい機会なので2ヶ月ほど Java にどっぷりハマってみました。 Apple の Swift に比べるとちょっと古い感じの言語ですが、似た部分が多くあり、 Swift がJava から色々取り入れているようで興味深いです。 もし、あなたがプログラミングの初心者ならば スッキリわかるJava入門 第2版 を1冊買って読みましょう。プログラミングをほとんど知らない人でも読めると思います。 私は2週間で読み終えました。Javaの基本とオブジェクト指向プログラミングの復習にもなって良かったです。関係者ではありませんが、宣伝しときます。 こっちはまだ読み終わってませんが、スッキリわかるJava入門 実践編 第2版 はもうちょっとレベルが上でラムダ式やデザインパターンにも触れています。Java についてネット上でも多くの解説がありますが、Java を使ったプログラミングをちゃんと学習するには、ちょっとお金を払って書籍を読むほうが近道のようです。

ベンチマークの結果

さて、前にベンチマークした100万回の行列の乗算を Java8 でも実行して Python、Lua、LuaJIT、C、アセンブラと速度を比較してみました。前のベンチマークの結果 と値が若干異なりますが、各言語で 4x4の正方行列同士の乗算を100万回行う時間の3回の平均値を同じ環境で取り直しました。

言語 備考 実行時間(秒)
Java 1.8.0 - 1.60
Python 2.7.3 リスト 154.06
NumPy 65.80
Python 3.2.3 リスト 157.78
NumPy 74.37
Lua 5.1.5 - 22.01
LuaJIT 2.0.0 テーブル(JIT) 2.41
テーブル(JIT OFF) 6.67
ffi double (JIT) 1.35
ffi double (JIT OFF) 72.24
gcc 4.6.3 -O0 2.28
-O2 0.66
アセンブリ ベクトル演算 0.64

Luajit で ffi を使った場合より少し遅くなりますが、Python より圧倒的に高速です。 Raspberry Pi でデスクトップ用のアプリをつくろうとすると Python + PyGTK と Java8 + Swing が有力ですが、速度を考えると Java8 + Swing ですね。 Android のアプリを作る練習にもなります。

java のコード

次のコードを「mulMatrix.java」というファイル名で保存して、コマンドラインから「javac mulMatrix.java」でコンパイルして mulMatrix.class を作り、その mulMatrix.class を 「java mulMatrix」で実行します。 java コマンドはクラス名を引数で与えますが、 クラスファイル名(mulMatrix.class) ではないことに注意して下さい。

  $ javac mulMatrix.java
  $ java mulMatrix
public class mulMatrix {

  public static void mulArrayMatrix(double[] mr, double[] ma, double[] mb) {
    mr[ 0] = ma[0] * mb[0] + ma[4] * mb[1] +  ma[8] * mb[2] + ma[12] * mb[3];
    mr[ 1] = ma[1] * mb[0] + ma[5] * mb[1] +  ma[9] * mb[2] + ma[13] * mb[3];
    mr[ 2] = ma[2] * mb[0] + ma[6] * mb[1] + ma[10] * mb[2] + ma[14] * mb[3];
    mr[ 3] = ma[3] * mb[0] + ma[7] * mb[1] + ma[11] * mb[2] + ma[15] * mb[3];
    mr[ 4] = ma[0] * mb[4] + ma[4] * mb[5] +  ma[8] * mb[6] + ma[12] * mb[7];
    mr[ 5] = ma[1] * mb[4] + ma[5] * mb[5] +  ma[9] * mb[6] + ma[13] * mb[7];
    mr[ 6] = ma[2] * mb[4] + ma[6] * mb[5] + ma[10] * mb[6] + ma[14] * mb[7];
    mr[ 7] = ma[3] * mb[4] + ma[7] * mb[5] + ma[11] * mb[6] + ma[15] * mb[7];
    mr[ 8] = ma[0] * mb[8] + ma[4] * mb[9] +  ma[8] * mb[10]+ ma[12] * mb[11];
    mr[ 9] = ma[1] * mb[8] + ma[5] * mb[9] +  ma[9] * mb[10]+ ma[13] * mb[11];
    mr[10] = ma[2] * mb[8] + ma[6] * mb[9] + ma[10] * mb[10]+ ma[14] * mb[11];
    mr[11] = ma[3] * mb[8] + ma[7] * mb[9] + ma[11] * mb[10]+ ma[15] * mb[11];
    mr[12] = ma[0] * mb[12]+ ma[4] * mb[13]+  ma[8] * mb[14]+ ma[12] * mb[15];
    mr[13] = ma[1] * mb[12]+ ma[5] * mb[13]+  ma[9] * mb[14]+ ma[13] * mb[15];
    mr[14] = ma[2] * mb[12]+ ma[6] * mb[13]+ ma[10] * mb[14]+ ma[14] * mb[15];
    mr[15] = ma[3] * mb[12]+ ma[7] * mb[13]+ ma[11] * mb[14]+ ma[15] * mb[15];
  }

  public static void printArrayMatrix4(double[] ma) {
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
                       ma[0], ma[4], ma[8],  ma[12]);
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
           ma[1], ma[5], ma[9],  ma[13]);
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
           ma[2], ma[6], ma[10], ma[14]);
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
           ma[3], ma[7], ma[11], ma[15]);
  }

  public static void main(String[] args) {
    double[] A = {1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4};
    double[] B = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
    double[] C = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    long start = System.nanoTime();

    for (int i=0; i<1000000; i++) {
      mulArrayMatrix(C, A, B);
    }

    long _end = System.nanoTime();
    System.out.printf("Elapsed Time = %10.8f sec\n",
                      (_end - start)/1000000000.0 );
    printArrayMatrix4(C);
  }
}

ベンチマークを実行したカーネルのバージョンと java のコンパイラのバージョンをここに記録しておきます。

$ javac -version
javac 1.8.0

$ uname -a
Linux raspberrypi 3.12.35+ #730 PREEMPT Fri Dec 19 18:31:24 GMT 2014 armv6l GNU/Linux

乗算方法の検討

Raspberry Pi上で java も Luajit と匹敵するほど速いことが分かりました。 4x4の正方行列の乗算の計算方法を工夫して、もう少し最適化できないか試してみました。

[1] Array element

言語間の速度比較で使ったコードです。配列からなる行列の要素をそのまま乗算に使っています。配列の各要素は4回使われるため、配列のインデックスによる参照が遅いと不利な方法です。

  public static void mulArrayMatrix(double[] mr, double[] ma, double[] mb)
  {
    mr[ 0] = ma[0] * mb[0] + ma[4] * mb[1] +  ma[8] * mb[2] + ma[12] * mb[3];
    mr[ 1] = ma[1] * mb[0] + ma[5] * mb[1] +  ma[9] * mb[2] + ma[13] * mb[3];
    mr[ 2] = ma[2] * mb[0] + ma[6] * mb[1] + ma[10] * mb[2] + ma[14] * mb[3];
    mr[ 3] = ma[3] * mb[0] + ma[7] * mb[1] + ma[11] * mb[2] + ma[15] * mb[3];
    mr[ 4] = ma[0] * mb[4] + ma[4] * mb[5] +  ma[8] * mb[6] + ma[12] * mb[7];
    mr[ 5] = ma[1] * mb[4] + ma[5] * mb[5] +  ma[9] * mb[6] + ma[13] * mb[7];
    mr[ 6] = ma[2] * mb[4] + ma[6] * mb[5] + ma[10] * mb[6] + ma[14] * mb[7];
    mr[ 7] = ma[3] * mb[4] + ma[7] * mb[5] + ma[11] * mb[6] + ma[15] * mb[7];
    mr[ 8] = ma[0] * mb[8] + ma[4] * mb[9] +  ma[8] * mb[10]+ ma[12] * mb[11];
    mr[ 9] = ma[1] * mb[8] + ma[5] * mb[9] +  ma[9] * mb[10]+ ma[13] * mb[11];
    mr[10] = ma[2] * mb[8] + ma[6] * mb[9] + ma[10] * mb[10]+ ma[14] * mb[11];
    mr[11] = ma[3] * mb[8] + ma[7] * mb[9] + ma[11] * mb[10]+ ma[15] * mb[11];
    mr[12] = ma[0] * mb[12]+ ma[4] * mb[13]+  ma[8] * mb[14]+ ma[12] * mb[15];
    mr[13] = ma[1] * mb[12]+ ma[5] * mb[13]+  ma[9] * mb[14]+ ma[13] * mb[15];
    mr[14] = ma[2] * mb[12]+ ma[6] * mb[13]+ ma[10] * mb[14]+ ma[14] * mb[15];
    mr[15] = ma[3] * mb[12]+ ma[7] * mb[13]+ ma[11] * mb[14]+ ma[15] * mb[15];
  }

[2] FOR Loop

行列の乗算を3重のループで行うコードです。短いコードですが、最も内側の命令「mr[i*4 + j]+= mb[i*4 + k] * ma[k*4 + j]」は 64回実行されます。

  public static void mulArrayMatrix2(double[] mr, double[] ma, double[] mb)
  {
    for(int i=0; i<4; i++)
      for(int j=0; j<4; j++) {
  mr[i*4 + j] = 0.0;
        for (int k=0; k<4; k++) 
          mr[i*4 + j]+= mb[i*4 + k] * ma[k*4 + j];
      }
  }

[3] Local variables

ローカル変数に配列の要素をコピーすることで、配列のインデックスによる参照を1回だけにしています。ローカル変数の割り当てと代入の速度が配列の参照より速ければ結果は良くなるはずです。

  public static void mulArrayMatrix3(double[] mr, double[] ma, double[] mb) {
    double ma0 = ma[0]; 
    double ma1 = ma[1]; 
    double ma2 = ma[2]; 
    double ma3 = ma[3]; 
    double ma4 = ma[4]; 
    double ma5 = ma[5]; 
    double ma6 = ma[6]; 
    double ma7 = ma[7]; 
    double ma8 = ma[8]; 
    double ma9 = ma[9]; 
    double ma10 = ma[10]; 
    double ma11 = ma[11]; 
    double ma12 = ma[12]; 
    double ma13 = ma[13]; 
    double ma14 = ma[14]; 
    double ma15 = ma[15]; 
    double mb0 = mb[0]; 
    double mb1 = mb[1]; 
    double mb2 = mb[2]; 
    double mb3 = mb[3]; 
    double mb4 = mb[4]; 
    double mb5 = mb[5]; 
    double mb6 = mb[6]; 
    double mb7 = mb[7]; 
    double mb8 = mb[8]; 
    double mb9 = mb[9]; 
    double mb10 = mb[10]; 
    double mb11 = mb[11]; 
    double mb12 = mb[12]; 
    double mb13 = mb[13]; 
    double mb14 = mb[14]; 
    double mb15 = mb[15]; 

    mr[ 0] = ma0 * mb0 + ma4 * mb1 +  ma8 * mb2 + ma12 * mb3;
    mr[ 1] = ma1 * mb0 + ma5 * mb1 +  ma9 * mb2 + ma13 * mb3;
    mr[ 2] = ma2 * mb0 + ma6 * mb1 + ma10 * mb2 + ma14 * mb3;
    mr[ 3] = ma3 * mb0 + ma7 * mb1 + ma11 * mb2 + ma15 * mb3;
    mr[ 4] = ma0 * mb4 + ma4 * mb5 +  ma8 * mb6 + ma12 * mb7;
    mr[ 5] = ma1 * mb4 + ma5 * mb5 +  ma9 * mb6 + ma13 * mb7;
    mr[ 6] = ma2 * mb4 + ma6 * mb5 + ma10 * mb6 + ma14 * mb7;
    mr[ 7] = ma3 * mb4 + ma7 * mb5 + ma11 * mb6 + ma15 * mb7;
    mr[ 8] = ma0 * mb8 + ma4 * mb9 +  ma8 * mb10+ ma12 * mb11;
    mr[ 9] = ma1 * mb8 + ma5 * mb9 +  ma9 * mb10+ ma13 * mb11;
    mr[10] = ma2 * mb8 + ma6 * mb9 + ma10 * mb10+ ma14 * mb11;
    mr[11] = ma3 * mb8 + ma7 * mb9 + ma11 * mb10+ ma15 * mb11;
    mr[12] = ma0 * mb12+ ma4 * mb13+  ma8 * mb14+ ma12 * mb15;
    mr[13] = ma1 * mb12+ ma5 * mb13+  ma9 * mb14+ ma13 * mb15;
    mr[14] = ma2 * mb12+ ma6 * mb13+ ma10 * mb14+ ma14 * mb15;
    mr[15] = ma3 * mb12+ ma7 * mb13+ ma11 * mb14+ ma15 * mb15;
  }

[4] Static properties

グローバル変数のように、クラスにスタティックなプロパティを用意しておいて、毎回ローカル変数を割り当てるのを防いでいます。普通はこんなコードは書きませんが、まあ比較ということで。

  static double ma0,ma1,ma2,ma3,ma4,ma5,ma6,ma7,ma8,ma9,ma10,
                ma11,ma12,ma13,ma14,ma15;
  static double mb0,mb1,mb2,mb3,mb4,mb5,mb6,mb7,mb8,mb9,mb10,
                mb11,mb12,mb13,mb14,mb15;
  public static void mulArrayMatrix4(double[] mr, double[] ma, double[] mb)
  {
    ma0 = ma[0]; 
    ma1 = ma[1]; 
    ma2 = ma[2]; 
    ma3 = ma[3]; 
    ma4 = ma[4]; 
    ma5 = ma[5]; 
    ma6 = ma[6]; 
    ma7 = ma[7]; 
    ma8 = ma[8]; 
    ma9 = ma[9]; 
    ma10 = ma[10]; 
    ma11 = ma[11]; 
    ma12 = ma[12]; 
    ma13 = ma[13]; 
    ma14 = ma[14]; 
    ma15 = ma[15]; 
    mb0 = mb[0]; 
    mb1 = mb[1]; 
    mb2 = mb[2]; 
    mb3 = mb[3]; 
    mb4 = mb[4]; 
    mb5 = mb[5]; 
    mb6 = mb[6]; 
    mb7 = mb[7]; 
    mb8 = mb[8]; 
    mb9 = mb[9]; 
    mb10 = mb[10]; 
    mb11 = mb[11]; 
    mb12 = mb[12]; 
    mb13 = mb[13]; 
    mb14 = mb[14]; 
    mb15 = mb[15]; 

    mr[ 0] = ma0 * mb0 + ma4 * mb1 +  ma8 * mb2 + ma12 * mb3;
    mr[ 1] = ma1 * mb0 + ma5 * mb1 +  ma9 * mb2 + ma13 * mb3;
    mr[ 2] = ma2 * mb0 + ma6 * mb1 + ma10 * mb2 + ma14 * mb3;
    mr[ 3] = ma3 * mb0 + ma7 * mb1 + ma11 * mb2 + ma15 * mb3;
    mr[ 4] = ma0 * mb4 + ma4 * mb5 +  ma8 * mb6 + ma12 * mb7;
    mr[ 5] = ma1 * mb4 + ma5 * mb5 +  ma9 * mb6 + ma13 * mb7;
    mr[ 6] = ma2 * mb4 + ma6 * mb5 + ma10 * mb6 + ma14 * mb7;
    mr[ 7] = ma3 * mb4 + ma7 * mb5 + ma11 * mb6 + ma15 * mb7;
    mr[ 8] = ma0 * mb8 + ma4 * mb9 +  ma8 * mb10+ ma12 * mb11;
    mr[ 9] = ma1 * mb8 + ma5 * mb9 +  ma9 * mb10+ ma13 * mb11;
    mr[10] = ma2 * mb8 + ma6 * mb9 + ma10 * mb10+ ma14 * mb11;
    mr[11] = ma3 * mb8 + ma7 * mb9 + ma11 * mb10+ ma15 * mb11;
    mr[12] = ma0 * mb12+ ma4 * mb13+  ma8 * mb14+ ma12 * mb15;
    mr[13] = ma1 * mb12+ ma5 * mb13+  ma9 * mb14+ ma13 * mb15;
    mr[14] = ma2 * mb12+ ma6 * mb13+ ma10 * mb14+ ma14 * mb15;
    mr[15] = ma3 * mb12+ ma7 * mb13+ ma11 * mb14+ ma15 * mb15;
  }

[5] java.nio.DoubleBuffer

[2] と同じように行列の乗算を3重のループで行うコードですが、配列の代わりに java.nio.DoubleBuffer を使ったものです。 luajit では ffi を使ってネイティブなバッファーを使った場合が最も速かったので、java でも同じような感じで試してみました。

  public static void mulArrayMatrix5(DoubleBuffer fr,
                                     DoubleBuffer fa, DoubleBuffer fb)
  {
    for(int i=0; i<4; i++)
      for(int j=0; j<4; j++) {
  fr.put(i*4 + j, 0.0);
        for (int k=0; k<4; k++)
          fr.put(i*4 + j, fr.get(i*4 + j) + fb.get(i*4 + k) * fa.get(k*4 + j));
      }
  }

結果のまとめ

次の表は、Raspberry Pi の Java8 の環境でコンパイルと実行をしたもの、Java7でコンパイルしたものを Raspberry Pi の Java8 の環境で実行したもの、MacOSX(i7 2.2GHz) の Java7 の環境でコンパイルと実行をしたものを比較しています。 Raspberry Pi(ARMv6 800MHz) の Java8 の環境でコンパイルと実行をしたものと、MacOSX(i7:x86_64 2.2GHz) の Java7 の環境でコンパイルと実行をしたものでは結果が大きく異なっています。 Java7 と Java8 で異なるのではなく、ARM と x86_64 の違いによるものでしょう。x86_64 では速度の差が小さいことから最適化が進んでいると思います。 Raspberry Pi ではループを使った方法ではかなり遅くなるので避けたほうが良いようです。

- java8c / java8r [ARM] java7c / java8r [ARM] java7c / java7r [i7]
[1] Array elements 1.626 秒 100 % 1.636 秒 100 % 0.055 秒 100 %
[2] FOR Loop 5.610 秒 345 % 5.629 秒 344 % 0.050 秒 90 %
[3] Local variables 1.719 秒 106 % 1.718 秒 105 % 0.045 秒 82 %
[4] Static properties 1.953 秒 120 % 1.958 秒 120 % 0.076 秒 139 %
[5] java.nio.DoubleBuffer 20.933 秒 1288 % 20.955 秒 1281 % 0.066 秒 120 %

ソースコード

// mulMatrix2.java
//   2014/12/28
//   Copyright (c) 2014 Jun Mizutani,
//   released under the MIT open source license.

import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.nio.DoubleBuffer;

public class mulMatrix2 {
  private static DoubleBuffer fb;
  private static DoubleBuffer fa;
  private static DoubleBuffer fr;

  // [1] Array element
  public static void mulArrayMatrix(double[] mr, double[] ma, double[] mb)
  {
    mr[ 0] = ma[0] * mb[0] + ma[4] * mb[1] +  ma[8] * mb[2] + ma[12] * mb[3];
    mr[ 1] = ma[1] * mb[0] + ma[5] * mb[1] +  ma[9] * mb[2] + ma[13] * mb[3];
    mr[ 2] = ma[2] * mb[0] + ma[6] * mb[1] + ma[10] * mb[2] + ma[14] * mb[3];
    mr[ 3] = ma[3] * mb[0] + ma[7] * mb[1] + ma[11] * mb[2] + ma[15] * mb[3];
    mr[ 4] = ma[0] * mb[4] + ma[4] * mb[5] +  ma[8] * mb[6] + ma[12] * mb[7];
    mr[ 5] = ma[1] * mb[4] + ma[5] * mb[5] +  ma[9] * mb[6] + ma[13] * mb[7];
    mr[ 6] = ma[2] * mb[4] + ma[6] * mb[5] + ma[10] * mb[6] + ma[14] * mb[7];
    mr[ 7] = ma[3] * mb[4] + ma[7] * mb[5] + ma[11] * mb[6] + ma[15] * mb[7];
    mr[ 8] = ma[0] * mb[8] + ma[4] * mb[9] +  ma[8] * mb[10]+ ma[12] * mb[11];
    mr[ 9] = ma[1] * mb[8] + ma[5] * mb[9] +  ma[9] * mb[10]+ ma[13] * mb[11];
    mr[10] = ma[2] * mb[8] + ma[6] * mb[9] + ma[10] * mb[10]+ ma[14] * mb[11];
    mr[11] = ma[3] * mb[8] + ma[7] * mb[9] + ma[11] * mb[10]+ ma[15] * mb[11];
    mr[12] = ma[0] * mb[12]+ ma[4] * mb[13]+  ma[8] * mb[14]+ ma[12] * mb[15];
    mr[13] = ma[1] * mb[12]+ ma[5] * mb[13]+  ma[9] * mb[14]+ ma[13] * mb[15];
    mr[14] = ma[2] * mb[12]+ ma[6] * mb[13]+ ma[10] * mb[14]+ ma[14] * mb[15];
    mr[15] = ma[3] * mb[12]+ ma[7] * mb[13]+ ma[11] * mb[14]+ ma[15] * mb[15];
  }

  // [2] FOR Loop
  public static void mulArrayMatrix2(double[] mr, double[] ma, double[] mb)
  {
    for(int i=0; i<4; i++)
      for(int j=0; j<4; j++) {
  mr[i*4 + j] = 0.0;
        for (int k=0; k<4; k++) 
          mr[i*4 + j]+= mb[i*4 + k] * ma[k*4 + j];
      }
  }

  // [3] Local variables
  public static void mulArrayMatrix3(double[] mr, double[] ma, double[] mb) {
    double ma0 = ma[0]; 
    double ma1 = ma[1]; 
    double ma2 = ma[2]; 
    double ma3 = ma[3]; 
    double ma4 = ma[4]; 
    double ma5 = ma[5]; 
    double ma6 = ma[6]; 
    double ma7 = ma[7]; 
    double ma8 = ma[8]; 
    double ma9 = ma[9]; 
    double ma10 = ma[10]; 
    double ma11 = ma[11]; 
    double ma12 = ma[12]; 
    double ma13 = ma[13]; 
    double ma14 = ma[14]; 
    double ma15 = ma[15]; 
    double mb0 = mb[0]; 
    double mb1 = mb[1]; 
    double mb2 = mb[2]; 
    double mb3 = mb[3]; 
    double mb4 = mb[4]; 
    double mb5 = mb[5]; 
    double mb6 = mb[6]; 
    double mb7 = mb[7]; 
    double mb8 = mb[8]; 
    double mb9 = mb[9]; 
    double mb10 = mb[10]; 
    double mb11 = mb[11]; 
    double mb12 = mb[12]; 
    double mb13 = mb[13]; 
    double mb14 = mb[14]; 
    double mb15 = mb[15]; 

    mr[ 0] = ma0 * mb0 + ma4 * mb1 +  ma8 * mb2 + ma12 * mb3;
    mr[ 1] = ma1 * mb0 + ma5 * mb1 +  ma9 * mb2 + ma13 * mb3;
    mr[ 2] = ma2 * mb0 + ma6 * mb1 + ma10 * mb2 + ma14 * mb3;
    mr[ 3] = ma3 * mb0 + ma7 * mb1 + ma11 * mb2 + ma15 * mb3;
    mr[ 4] = ma0 * mb4 + ma4 * mb5 +  ma8 * mb6 + ma12 * mb7;
    mr[ 5] = ma1 * mb4 + ma5 * mb5 +  ma9 * mb6 + ma13 * mb7;
    mr[ 6] = ma2 * mb4 + ma6 * mb5 + ma10 * mb6 + ma14 * mb7;
    mr[ 7] = ma3 * mb4 + ma7 * mb5 + ma11 * mb6 + ma15 * mb7;
    mr[ 8] = ma0 * mb8 + ma4 * mb9 +  ma8 * mb10+ ma12 * mb11;
    mr[ 9] = ma1 * mb8 + ma5 * mb9 +  ma9 * mb10+ ma13 * mb11;
    mr[10] = ma2 * mb8 + ma6 * mb9 + ma10 * mb10+ ma14 * mb11;
    mr[11] = ma3 * mb8 + ma7 * mb9 + ma11 * mb10+ ma15 * mb11;
    mr[12] = ma0 * mb12+ ma4 * mb13+  ma8 * mb14+ ma12 * mb15;
    mr[13] = ma1 * mb12+ ma5 * mb13+  ma9 * mb14+ ma13 * mb15;
    mr[14] = ma2 * mb12+ ma6 * mb13+ ma10 * mb14+ ma14 * mb15;
    mr[15] = ma3 * mb12+ ma7 * mb13+ ma11 * mb14+ ma15 * mb15;
  }

  // [4] Static properties
  static double ma0,ma1,ma2,ma3,ma4,ma5,ma6,ma7,ma8,ma9,ma10,
                ma11,ma12,ma13,ma14,ma15;
  static double mb0,mb1,mb2,mb3,mb4,mb5,mb6,mb7,mb8,mb9,mb10,
                mb11,mb12,mb13,mb14,mb15;
  public static void mulArrayMatrix4(double[] mr, double[] ma, double[] mb)
  {
    ma0 = ma[0]; 
    ma1 = ma[1]; 
    ma2 = ma[2]; 
    ma3 = ma[3]; 
    ma4 = ma[4]; 
    ma5 = ma[5]; 
    ma6 = ma[6]; 
    ma7 = ma[7]; 
    ma8 = ma[8]; 
    ma9 = ma[9]; 
    ma10 = ma[10]; 
    ma11 = ma[11]; 
    ma12 = ma[12]; 
    ma13 = ma[13]; 
    ma14 = ma[14]; 
    ma15 = ma[15]; 
    mb0 = mb[0]; 
    mb1 = mb[1]; 
    mb2 = mb[2]; 
    mb3 = mb[3]; 
    mb4 = mb[4]; 
    mb5 = mb[5]; 
    mb6 = mb[6]; 
    mb7 = mb[7]; 
    mb8 = mb[8]; 
    mb9 = mb[9]; 
    mb10 = mb[10]; 
    mb11 = mb[11]; 
    mb12 = mb[12]; 
    mb13 = mb[13]; 
    mb14 = mb[14]; 
    mb15 = mb[15]; 

    mr[ 0] = ma0 * mb0 + ma4 * mb1 +  ma8 * mb2 + ma12 * mb3;
    mr[ 1] = ma1 * mb0 + ma5 * mb1 +  ma9 * mb2 + ma13 * mb3;
    mr[ 2] = ma2 * mb0 + ma6 * mb1 + ma10 * mb2 + ma14 * mb3;
    mr[ 3] = ma3 * mb0 + ma7 * mb1 + ma11 * mb2 + ma15 * mb3;
    mr[ 4] = ma0 * mb4 + ma4 * mb5 +  ma8 * mb6 + ma12 * mb7;
    mr[ 5] = ma1 * mb4 + ma5 * mb5 +  ma9 * mb6 + ma13 * mb7;
    mr[ 6] = ma2 * mb4 + ma6 * mb5 + ma10 * mb6 + ma14 * mb7;
    mr[ 7] = ma3 * mb4 + ma7 * mb5 + ma11 * mb6 + ma15 * mb7;
    mr[ 8] = ma0 * mb8 + ma4 * mb9 +  ma8 * mb10+ ma12 * mb11;
    mr[ 9] = ma1 * mb8 + ma5 * mb9 +  ma9 * mb10+ ma13 * mb11;
    mr[10] = ma2 * mb8 + ma6 * mb9 + ma10 * mb10+ ma14 * mb11;
    mr[11] = ma3 * mb8 + ma7 * mb9 + ma11 * mb10+ ma15 * mb11;
    mr[12] = ma0 * mb12+ ma4 * mb13+  ma8 * mb14+ ma12 * mb15;
    mr[13] = ma1 * mb12+ ma5 * mb13+  ma9 * mb14+ ma13 * mb15;
    mr[14] = ma2 * mb12+ ma6 * mb13+ ma10 * mb14+ ma14 * mb15;
    mr[15] = ma3 * mb12+ ma7 * mb13+ ma11 * mb14+ ma15 * mb15;
  }

  // [5] java.nio.DoubleBuffer
  public static void mulArrayMatrix5(DoubleBuffer fr,
                                     DoubleBuffer fa, DoubleBuffer fb)
  {
    for(int i=0; i<4; i++)
      for(int j=0; j<4; j++) {
  fr.put(i*4 + j, 0.0);
        for (int k=0; k<4; k++)
          fr.put(i*4 + j, fr.get(i*4 + j) + fb.get(i*4 + k) * fa.get(k*4 + j));
      }
  }

  public static void printArrayMatrix4(double[] ma) {
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
                       ma[0], ma[4], ma[8],  ma[12]);
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
           ma[1], ma[5], ma[9],  ma[13]);
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
           ma[2], ma[6], ma[10], ma[14]);
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
           ma[3], ma[7], ma[11], ma[15]);
  }

   public static void printBufferMatrix4(DoubleBuffer d) {
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
                       d.get(0), d.get(4), d.get(8),  d.get(12));
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
                       d.get(1), d.get(5), d.get(9),  d.get(13));
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
                       d.get(2), d.get(6), d.get(10),  d.get(14));
    System.out.printf(" %10.5e  %10.5e  %10.5e  %10.5e\n",
                       d.get(3), d.get(7), d.get(11),  d.get(15));
  }

  public static void printTime(long start)
  {
    long _end = System.nanoTime();
    _end = System.nanoTime();
    System.out.printf("Elapsed Time = %10.8f sec\n",
                      (_end - start)/1000000000.0 );
  }

  public static void main(String[] args) {
    double[] A = {1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4};
    double[] B = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
    double[] C = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    ByteBuffer bb;
    
    // Initialize DoubleBuffers
    bb = ByteBuffer.allocateDirect(16 * 8).order(ByteOrder.nativeOrder());
    fa = bb.asDoubleBuffer();
    bb = ByteBuffer.allocateDirect(16 * 8).order(ByteOrder.nativeOrder());
    fb = bb.asDoubleBuffer();
    bb = ByteBuffer.allocateDirect(16 * 8).order(ByteOrder.nativeOrder());
    fr = bb.asDoubleBuffer();
    for (int i=0; i<16; i++) {
      fa.put(i,  A[i]);
      fb.put(i,  B[i]);
      fr.put(i,  C[i]);
    }

    // [1] Array elements
    System.out.printf("[1] Array elements\n");
    long start = System.nanoTime();
    for (int i=0; i<1000000; i++) {
      mulArrayMatrix(C, A, B);
    }
    printTime(start);
    printArrayMatrix4(C);

    // [2] FOR Loop
    System.out.printf("[2] FOR Loop\n");
    start = System.nanoTime();
    for (int i=0; i<1000000; i++) {
      mulArrayMatrix2(C, A, B);
    }
    printTime(start);
    printArrayMatrix4(C);          

    // [3] Local variables
    System.out.printf("[3] Local variables\n");
    start = System.nanoTime();
    for (int i=0; i<1000000; i++) {
      mulArrayMatrix3(C, A, B);
    }
    printTime(start);
    printArrayMatrix4(C);

    // [4] Static properties
    System.out.printf("[4] Static properties\n");
    start = System.nanoTime();
    for (int i=0; i<1000000; i++) {
      mulArrayMatrix4(C, A, B);
    }
    printTime(start);
    printArrayMatrix4(C);          

    // [5] java.nio.DoubleBuffer
    System.out.printf("[5] java.nio.DoubleBuffer\n");
    start = System.nanoTime();
    for (int i=0; i<1000000; i++) {
      mulArrayMatrix5(fr, fa, fb);
    }
    printTime(start);
    printBufferMatrix4(fr);          
  }
}
$ javac mulMatrix2.java
$ java mulMatrix2
[1] Array elements
Elapsed Time = 1.62580848 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
[2] FOR Loop
Elapsed Time = 5.61011917 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
[3] Local variables
Elapsed Time = 1.71906921 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
[4] Static properties
Elapsed Time = 1.95252907 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
[5] java.nio.DoubleBuffer
Elapsed Time = 20.93300011 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

さらなる高速化

4x4の正方行列を3次元コンピュータグラフィックで使う場合は、通常、4行目が [0.0, 0.0, 0.0, 1.0] となる行列(m[3]、m[7]、m[11]が0.0、m[15]が1.0) であることを利用すると、次のように手抜きできます。

  public static void mulArrayMatrix(double[] mr, double[] ma, double[] mb) {
    mr[ 0] = ma[0] * mb[0] + ma[4] * mb[1] +  ma[8] * mb[2] ;
    mr[ 1] = ma[1] * mb[0] + ma[5] * mb[1] +  ma[9] * mb[2] ;
    mr[ 2] = ma[2] * mb[0] + ma[6] * mb[1] + ma[10] * mb[2] ;
    mr[ 3] = 0.0;
    mr[ 4] = ma[0] * mb[4] + ma[4] * mb[5] +  ma[8] * mb[6] ;
    mr[ 5] = ma[1] * mb[4] + ma[5] * mb[5] +  ma[9] * mb[6] ;
    mr[ 6] = ma[2] * mb[4] + ma[6] * mb[5] + ma[10] * mb[6] ;
    mr[ 7] = 0.0;
    mr[ 8] = ma[0] * mb[8] + ma[4] * mb[9] +  ma[8] * mb[10];
    mr[ 9] = ma[1] * mb[8] + ma[5] * mb[9] +  ma[9] * mb[10];
    mr[10] = ma[2] * mb[8] + ma[6] * mb[9] + ma[10] * mb[10];
    mr[11] = 0.0;
    mr[12] = ma[0] * mb[12]+ ma[4] * mb[13]+  ma[8] * mb[14];
    mr[13] = ma[1] * mb[12]+ ma[5] * mb[13]+  ma[9] * mb[14];
    mr[14] = ma[2] * mb[12]+ ma[6] * mb[13]+ ma[10] * mb[14];
    mr[15] = 1.0;
  }

実行すると:

Elapsed Time = 1.00589521 sec
 3.00000e+00  1.50000e+01  2.70000e+01  3.90000e+01
 6.00000e+00  3.00000e+01  5.40000e+01  7.80000e+01
 9.00000e+00  4.50000e+01  8.10000e+01  1.17000e+02
 0.00000e+00  0.00000e+00  0.00000e+00  1.00000e+00

倍精度の4x4行列の乗算 100万回が 1.0秒で終わりました。十分な速度です。



「Raspberry Pi でラムダ式」 に続く...