rust ndarrayの行列演算をマルチスレッドで行いたい

Question

やりたいこと

rust ndarrayを用いて線形代数や一般的な行列演算を高速に行いたいと思っています。高速に行いたいため、効率よくCPUのコアを使って欲しいと思っています。

やったこと

ndarrayとndarray_linalgを用いて、行列演算を試験的逆行列を求めるコードかき、topコマンドでCPUの使用率をモニターしました。
使用したコードは以下です。

use ndarray::{Array2, Array};
use ndarray_linalg::InverseInto;
use ndarray_rand::RandomExt;
use ndarray_rand::rand_distr::Normal;

fn main() {
    let a = Array::random((100000,100000), Normal::new(1.,1.).unwrap());
    let inv_a = a.inv_into().unwrap();
}

結果

CPUの使用率をモニターしたところ100%に張り付いており1コアのみで演算を行っているようにおもいました。

調べたこと

ndarray_linalgのソースを読んで調べてみたのですが、逆行列を求めている該当の部分に関しては、ndarray-linalg/src/solve.rsの部分だと考えています。

impl<A, S> InverseInto for ArrayBase<S, Ix2>
where
    A: Scalar + Lapack,
    S: DataMut<Elem = A> + RawDataClone,
{
    type Output = Self;

    fn inv_into(self) -> Result<Self::Output> {
        let f = self.factorize_into()?;
        f.inv_into()
    }
}

とあり、ScalarかLapackにinv_intoが実装されていると考えました。そこで、Scalarはcauchyで実装されているもので、inv_intoが実装されていなさそうだとおもいました。
そのため、Lapackで逆行列を求めるinv_intoが実装されていると考えました。
Lapackは/ndarray-linalg/lax/src/lib.rsに実装されておりました。そのため、inv_intoはOperatorNorm_,
QR_,SVD_,SVDDC_,Solve_,Solveh_,Cholesky_,Eig_,Eigh_,Triangular_,Tridiagonal_,Rcond_,LeastSquaresSvdDivideConquer_のうちのどれかに実装されていそうだと考えました。結局、Solve_で実装されていると考えました。invの実装しか見つからなかったのでinv_intoの中身がどこにあるのかは結局よくわかりませんでした。

また、invの中身に関しても、Rust初心者なので理解できませんでした。
結局、なぜマルチスレッドで計算できないのかがよくわかりませんでした。

お答えいただきたいこと

intel_mklを用いて計算しているのであればマルチスレッドで計算されると考えておりました。
なぜマルチスレッドにならないか、どのようにすれば解決するのかを教えていただきたいです。

Answer 1

まだリリースしていない 0.13.0 (GitHub の master) ベースでの話になってしまって恐縮ですが,
ndarray-linalg から Intel MKL を呼び出す場合 intel-mkl feature を指定しますが、これは MKL の serial 実装 mkl-static-lp64-seq をリンクするため並列化されていません。MKL自体に並列実行させるためには mkl-static-lp64-iomp をリンクする必要がありますが、linalg 側にこれを指定するオプションは現在作っていないので Cargo.toml の [patch] で intel-mkl-src crate 自体を置き換えるか、linalg 自体を変更するなら lax/Cargo.toml ここの seq を iomp に変更すると動くはずです

Comments

[BoKuToTuZenU]

ご回答いただきありがとうございます。こちらをご教授いただいた方法でやってみようかとおもいます。
また、今後のリリースでintel-mklの設定を変更できるようなfeatureを加える予定はありますか？

Answer 2

さわったことないからわからないですけど、
OMP_NUM_THREADS に値設定したらふえたりしないですかね？

https://jp.xlsoft.com/documents/intel/mkl/mac/Getting_Started_mac.htm

スレッド化
インテル MKL は、多くの個所でスレッド化されています。また、OpenMP* ソフトウェアを使用しています。
OpenMP ソフトウェアは、使用するスレッド数を設定する OMP_NUM_THREADS 環境変数を基にスレッド化を実行します。スレッド数は、プログラムを実行するシェルで設定することができます。
変数 OMP_NUM_THREADS が設定されていない場合、インテル MKL ソフトウェアはスレッド数を 1 として実行します。OMP_NUM_THREADS は、常にアプリケーションで使用するプロセッサー数に設定しておくことを推奨します。
1D FFT における 2 の累乗データについては、インテル MKL はインテル® Itanium® プロセッサー用にのみ並列化を行うことに注意してください。
スレッド数を変更するには、プログラムを実行するコマンドライン・シェルで次のように入力します。

export OMP_NUM_THREADS=<使用するスレッド数>

Comments

[BoKuToTuZenU]

ご回答いただき、ありがとうございます。export OMP_NUM_THREADS=12と設定して計算させてみましたが、1スレッドで計算を行っているようです。また、私がリポジトリにあるコードを読んだ感じではどこでintel mklやOpenblasなどのライブラリを読み出しているのかがわかりませんでした。