NetBeans IDE で JLAPACK を試す (2)

JLAPACK ライブラリが NetBeans IDE 上で利用できるようになったので、以前、gfortran + LAPACK で試したサンプルを、Java + JLAPACK でも試してみました。参考サイト a. で、LAPACK の最小二乗問題を扱う倍精度用の手続き副プログラム DGELS ^b. を利用した QR 分解のサンプルを紹介していますので、これを Java 用に書き換えました。

サンプルは、100000 組の x と y が対になっているデータ fdata01.zip (fdata01.txt) を読み込んで、下記の4次の多項式に最小二乗法でフィッティングしたときの係数 c₁ 〜 c₅ を求めるというプログラムです。

y = c₀ + c₁x + c₂x² + c₃x³ + c₄x⁴

上記 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. というクラスも用意されているはずなのですが、このクラスがなぜか認識されなかったので、試すことができませんでした。

参考サイト

1. LAPACK の利用例 - Workshop Complex at bitWalk
2. Dgels - The Innovative Computing Laboratory (ICL) のサイトに掲載されている JLAPACK の文書
3. DGELS - 同上
4. NetBeans IDE で JLAPACK を試す - NetBeans IDE 上でライブラリの使い方を説明
5. bitWalk's workshop - JLAPACK