BlasSharp、という新しいライブラリを公開しました！

https://github.com/nuskey8/BlasSharp

GitHubの表示上はRustになっていますが、これはバインディングの生成に使っているだけで、本体は名前の通りC#向けのライブラリです。対象プラットフォームはパッケージにもよりますが、OpenBLASを使えば大体のWindows、macOS、Linux上で動作するはずです。

一応v0.2.0ということになっていますが、バインディング部分は既に安定しているため、十分使える状態にはなっているでしょう。

BLAS/LAPACK

そもそもBLAS/LAPACKとは何かというと、BLASは「Basic Linear Algebra Subprograms」、LAPACKは「Linear Algebra PACKage」の略称で、要するに線形代数計算のためのライブラリです。BLASは行列・ベクトルの基本的な操作(加算、スカラー倍、内積、行列積など)を提供するライブラリで、LAPACKはBLASを用いた連立一次方程式、固有値、固有ベクトルなどの計算を提供するライブラリになっています。

BLAS/LAPACKは元々Fortranで実装されたライブラリでしたが、今ではGotoBLASやその後継であるOpenBLAS、Intelが実装しているMKLなど、同じインターフェースを持つ様々な実装が公開されています。

行列の基礎的な演算程度なら自力で実装できそうに思えますが、基本的にはこれらのBLAS/LAPACKライブラリを利用することを強く推奨します。実際、NumPyなどの多くの数値計算ライブラリは上のようなBLAS/LAPACKのライブラリを組み込んで動作しています。

というのも線形代数計算は一般的に計算量が多く、普通に実装しただけでは全然速度が出ないんですね。例として、行列積(GEMM)を求める処理でベンチマークをとってみましょう。

using System;
using BenchmarkDotNet;
using BenchmarkDotNet.Attributes;
using BlasSharp;

namespace Benchmark;

public unsafe class DgemmBenchmark
{
    const int Size = 1024;

    double[] A, B, C;

    [GlobalSetup]
    public void Setup()
    {
        A = new double[Size * Size];
        B = new double[Size * Size];
        C = new double[Size * Size];

        var rand = new Random(0x1234);
        for (int i = 0; i < A.Length; i++) A[i] = rand.NextDouble();
        for (int i = 0; i < B.Length; i++) B[i] = rand.NextDouble();
    }

    [Benchmark]
    public void ForLoop()
    {
        // forループで普通に計算する
        for (int i = 0; i < Size; i++)
        {
            for (int j = 0; j < Size; j++)
            {
                double sum = 0;
                for (int k = 0; k < Size; k++)
                {
                    sum += A[i * Size + k] * B[k * Size + j];
                }
                C[i * Size + j] = sum;
            }
        }
    }

    [Benchmark]
    public void OpenBlas()
    {
        // OpenBLASの実装を用いて計算する
        fixed (double* aPtr = A, bPtr = B, cPtr = C)
        {
            BlasSharp.OpenBlas.NativeMethods.cblas_dgemm(
                CBLAS_ORDER.CblasRowMajor, CBLAS_TRANSPOSE.CblasNoTrans, CBLAS_TRANSPOSE.CblasNoTrans,
                Size, Size, Size,
                1.0, aPtr, Size,
                bPtr, Size,
                0.0, cPtr, Size
            );
        }
    }

    [Benchmark]
    public void AppleAccelerate()
    {
        // AppleのAccelerateパッケージを用いて計算する
        fixed (double* aPtr = A, bPtr = B, cPtr = C)
        {
            BlasSharp.AppleAccelerate.NativeMethods.cblas_dgemm(
                CBLAS_ORDER.CblasRowMajor, CBLAS_TRANSPOSE.CblasNoTrans, CBLAS_TRANSPOSE.CblasNoTrans,
                Size, Size, Size,
                1.0, aPtr, Size,
                bPtr, Size,
                0.0, cPtr, Size
            );
        }
    }

}

手元のMacBook Pro(M2)で実行した結果は以下のようになります。

Method	Mean	Error	StdDev
ForLoop	1,423.217 ms	14.6145 ms	12.2038 ms
OpenBlas	16.703 ms	0.4721 ms	1.3546 ms
AppleAccelerate	6.781 ms	0.0585 ms	0.0457 ms

見ての通り、x85〜x200ほどの圧倒的な差が出ます。OpenBLASやAccelerateの実装はCPUレベルで限界までチューニングされており、普通の実装でこの速度を追い越すことは不可能です。そのため、通常はバインディングを介してこれらのライブラリを利用することになります。

BlasSharpの使い方

というわけで、BlasSharpではこれらのライブラリに対するバインディングを提供します。OpenBLAS、MKL、Accelerate(Apple)の3つのライブラリに対応しており、それぞれ別のパッケージとして提供しています。

using BlasSharp;
using static BlasSharp.OpenBlas.NativeMethods;

unsafe
{
    double[] A = [
        1.0, 2.0, 3.0,
        4.0, 5.0, 6.0,
        7.0, 8.0, 9.0
    ];

    double[] B = [
        1.0, 2.0, 3.0,
        4.0, 5.0, 6.0,
        7.0, 8.0, 9.0,
    ];

    double[] C = new double[9];

    fixed (double* a = A, b = B, c = C)
    {
        cblas_dgemm(CBLAS_ORDER.CblasColMajor,
            CBLAS_TRANSPOSE.CblasNoTrans,
            CBLAS_TRANSPOSE.CblasNoTrans,
            3, 3, 3, 1.0, a, 3, b, 3, 0.0, c, 3);
    }

    Console.WriteLine(string.Join(',', C));
}

もちろん、実行環境によってどのBLAS実装を選択すれば良いかは異なるでしょう。そこでBlasSharp.Sharedパッケージでは、BLASのAPIを抽象化するためのIBlasOperationsインターフェースを提供しています。BlasSharp.Sharedに含まれるのはインターフェースのみで、実際のIBlasOperationsの実装は上記のパッケージにはそれぞれのライブラリに対応したものが含まれています。

public interface IBlasOperations : IBlasLevel1, IBlasLevel2, IBlasLevel3
{
}

public unsafe interface IBlasLevel1
{
    void Saxpy(int n, float alpha, float* x, int incX, float* y, int incY);
    void Daxpy(int n, double alpha, double* x, int incX, double* y, int incY);
    void Caxpy(int n, void* alpha, void* x, int incX, void* y, int incY);
    void Zaxpy(int n, void* alpha, void* x, int incX, void* y, int incY);

    ...
}

BlasSharpは数値計算の基盤としての用途を想定しているため、特に高レベルAPIなどは用意していません。基本的にはunsafeでポインタを介して操作することになります。

とはいえSpan<T>で扱えた方が便利なので、拡張メソッドは用意しておこうかな...

rust-bindgen / csbindgenを用いたバインディング生成

BLAS/LAPACKの提供するルーチンは膨大な数になるため、FFIの部分を手書きするのは非現実的です。そのため何らかの方法でバインディングを自動生成することになるのですが、今回はrust-bindgenとcsbindgenを用いてRust経由で生成を行う方針にしています。

以下はIntelのMKLのバインディングを生成するbuild.rsです。

fn main() {
    bindgen::Builder::default()
        .headers(["headers/mkl_cblas.h", "headers/mkl_lapack.h"])
        .generate()
        .unwrap()
        .write_to_file("src/mkl.rs")
        .unwrap();

    csbindgen::Builder::default()
        .input_bindgen_file("src/mkl.rs")
        .rust_method_prefix("ffi_")
        .csharp_entry_point_prefix("")
        .csharp_method_prefix("")
        .csharp_disable_emit_dll_name(true)
        .csharp_namespace("BlasSharp.Mkl")
        .csharp_class_accessibility("public")
        .rust_file_header("use super::mkl::*;")
        .method_filter(|x| method_filter(x))
        .generate_to_file(
            "src/mkl_ffi.rs",
            "../../src/BlasSharp.Mkl/NativeMethods.g.cs",
        )
        .unwrap();
}

fn method_filter(x: String) -> bool {
    if x.starts_with("cblas_") {
        return true;
    };
    if x.starts_with("mkl_") {
        return true;
    };
    if x.ends_with("_") && !x.chars().any(|x| x.is_uppercase()) {
        return true;
    };
    return false;
}

bindgenがCのヘッダファイルを解析してRustコードを出力し、そのRustコードをcsbindgenが解析することでC#のコードを出力しています。csbindgenの使い方に関してはneueccさんの記事で詳しく解説がされているので、こちらを参考にすると良いでしょう。

実際めちゃくちゃ簡単にFFIが行えるのでおすすめです。csbindgenはいいぞ。

まとめ

そもそもC#で数値計算をやる人はおそらくほぼいない上、ほとんどただのバインディングでしかないので、正直これを直接使いたい人はあんまりいないんじゃないかな...という感じだと思います。実際私もより高度な数値計算ライブラリの開発を進めていて、BlasSharpはその基盤部分として利用するために開発したライブラリだったり。

とはいえ基盤として安定したものを提供できているとは思うので、もし使う機会があれば是非、ということで。

Discussion