JLAPACK ライブラリが NetBeans IDE 上で利用できるようになったので、以前、gfortran + LAPACK で試したサンプルを、Java + JLAPACK でも試してみました。参考サイト a. で、LAPACK の最小二乗問題を扱う倍精度用の手続き副プログラム DGELS b. を利用した QR 分解のサンプルを紹介していますので、これを Java 用に書き換えました。
サンプルは、100000 組の x と y が対になっているデータ fdata01.zip (fdata01.txt) を読み込んで、下記の4次の多項式に最小二乗法でフィッティングしたときの係数 c1 〜 c5 を求めるというプログラムです。
y = c0 + c1x + c2x2 + c3x3 + c4x4
上記 gfortran のサンプルプログラム ex101_01.f を Java に書き直した DgelsTest.java を以下に示します。
package dgelstest; import java.io.BufferedReader; import java.io.File; import java.io.FileNotFoundException; import java.io.FileReader; import java.io.IOException; import org.netlib.lapack.Dgels; import org.netlib.util.intW; public class DgelsTest { public static void main(String[] args) { intW info = new intW(2); int m = 100000, n = 5; int nrhs = 1, lda = m, ldb = m, lwork = 2 * n; double a[], b[], work[]; a = new double[m * n]; b = new double[m]; work = new double[lwork]; File file = new File("/home/bitwalk/work/ex101_01/fdata01.txt"); try (BufferedReader br = new BufferedReader(new FileReader(file))) { String str; int i = 0; while ((str = br.readLine()) != null) { double x = Double.parseDouble(str.substring(0, 6)); b[i] = Double.parseDouble(str.substring(6, 20)); for (int j = 0; j < n; j++) { a[i + m * j] = Math.pow(x, (double) j); } i++; } br.close(); } catch (FileNotFoundException e) { System.out.println(e); } catch (IOException e) { System.out.println(e); } Dgels.dgels("N", m, n, nrhs, a, 0, lda, b, 0, ldb, work, 0, lwork, info); if (info.val == 0) { System.out.println("# of DATA\t" + m); System.out.println("Coefficients"); for (int j = 0; j < n; j++) { System.out.println("C(" + j + ") = " + b[j]); } } else { System.out.println("Error!"); System.out.println("info = " + info.val); } } }
実行結果は以下のようになりました。
run: # of DATA 100000 Coefficients C(0) = 37.76585776132993 C(1) = 0.04968384703750248 C(2) = 6.988983959858456E-6 C(3) = -4.2174249160048426E-7 C(4) = 2.686562674530099E-10 ビルド成功(合計時間: 0秒)
いまどきの PC(と言っても、自分の PC は結構古いのですが)だと、Java を使っても、 100000 x 5 の行列計算程度では、データの読み込み時間を含めて大した時間がかかりません。
Dgels クラスでは、行列の要素を一次元配列に配置するので、C や Java の二次元配列における行と列の関係が、Fortran では逆になることを意識しなければなりません。なお、JLAPACK では、行列を二次元配列で扱える、DGELS c. というクラスも用意されているはずなのですが、このクラスがなぜか認識されなかったので、試すことができませんでした。
参考サイト
- LAPACK の利用例 - Workshop Complex at bitWalk
- Dgels - The Innovative Computing Laboratory (ICL) のサイトに掲載されている JLAPACK の文書
- DGELS - 同上
- NetBeans IDE で JLAPACK を試す - NetBeans IDE 上でライブラリの使い方を説明
- bitWalk's workshop - JLAPACK
0 件のコメント:
コメントを投稿