faerとnalgebraの比較

はじめに

これまでRustで有限要素法のプログラムを組んできましたが、行列の計算にはnalgebraの直接法を使ってきました。
ただ、nalgebraは疎行列の計算用ライブラリではないので計算速度がネックでした。
そんな中、faerというRustライブラリがあると知りました。疎行列の計算も入っているとのことでfaerとnalgebraで試計算用の行列を作って計算速度の比較を行いました。

使用したライブラリ

使用したライブラリは以下の通りです。

faer = { version = "0.24", features = ["rayon"] }
nalgebra = "0.35.0"
rayon = "1.12.0"

faerについてはfeatures = ["rayon"]と入れることで自動的に並列計算ができるようになるそうです。

行列とベクトルの作成

今回は有限要素法のような連立方程式を想定し、行列とベクトルの計算を行います。
有限要素法は一般に疎行列になるので疎行列と適当なベクトルで計算速度の比較を行いたいと思います。
疎行列は帯行列、ベクトルの中身は問わないとして以下のように作成しました。

let mut matrix_na = na::DMatrix::<f64>::zeros(N, N);
let mut matrix_fa = fa::Mat::<f64>::zeros(N, N);
let mut vector_na = na::DVector::<f64>::zeros(N);
let mut vector_fa = fa::Col::<f64>::zeros(N);
let bandwidth = 5;

for i in 0..N {
    for j in i.saturating_sub(bandwidth)..=(i + bandwidth).min(N - 1) {
        let v = if i == j {
            10.0
        } else {
            1.0 / ((i as isize - j as isize).abs() as f64 + 1.0)
        };

        matrix_na[(i, j)] = v;
        matrix_fa[(i, j)] = v;
    }

        let x = (i + 1) as f64;
        vector_na[i] = x;
        vector_fa[i] = x;
}

行列の計算

連立方程式の解法は直接法、反復法など数多く提案されていますが、ここでは一般的なLU分解を使いたいと思います。
ただし、faerでは疎行列の計算のため圧縮格納することにします。
faerの疎行列格納は一般的にCSC形式で保存されるそうで、ここでは一旦疎行列を作成し、それをCSC形式に圧縮して計算したいと思います。
なお、nalgebraは作成した疎行列を直接LU分解で計算します。
各ライブラリの計算コードは以下のようになります。

nalgebra

let _test_na = matrix_na
    .lu()
    .solve(&vector_na)
    .expect("Failed to solve linear equation");

faer

let (m, n) = matrix_fa.shape();
let mut triplets = Vec::<Triplet<usize, usize, f64>>::new();

for i in 0..m {
    for j in 0..n {
        let v = matrix_fa[(i, j)];

        if v != 0.0 {
            triplets.push(Triplet::new(i, j, v));
        }
    }
}

let sp_mat = SparseColMat::try_new_from_triplets(m, n, &triplets).unwrap();
let lu = sp_mat.sp_lu().unwrap();
let _test_fa = lu.solve(&vector_fa);

#　結果の比較
行列の大きさをいくつか変えて計算時間を比較してみました。
なお計算は全てM2 Macbookで行っています。
また、コンパイルの際はcargo run --releaseで最適化しています。

N=100

nalgebra: 12.797917ms
faer: 14.727292ms

N=1000

nalgebra: 121.752083ms
faer: 14.530875ms

N=10000

nalgebra: 96.210403s
faer: 284.440917ms

Nが増えると疎行列計算の効果がはっきり出ています。
また、faerだけを使ってN=100000の計算を行い並列処理ができたか確認しました。

RAYON_NUM_THREADS=1 cargo run --release
faer: 97.996636041s
RAYON_NUM_THREADS=8 cargo run --release
faer: 94.703522292s

計算結果はあまり変わらず並列処理ができていない印象です。
このあたりは詳しい方がいたらぜひご教示いただきたいです。

最後に

faerで連立方程式の計算が高速でできることが確認できました。
これを開発した有限要素解析プログラムに組み込みたいと思います。