Juliaの線形代数を超高速に!LinearAlgebra.BLAS徹底解説
Julia の LinearAlgebra
モジュールは、線形代数の基本的な操作(行列の積、行列式、逆行列、固有値計算など)を提供します。その中でも、LinearAlgebra.BLAS
は、BLAS (Basic Linear Algebra Subprograms) と呼ばれる、低レベルの線形代数演算の標準ライブラリへのインターフェースを提供します。
BLASは、ベクトルとベクトルの演算、行列とベクトルの演算、行列と行列の演算といった、線形代数で頻繁に登場する基本的な操作を、非常に高速に実行するために最適化されたルーチンの集合体です。これらのルーチンは通常、アセンブリ言語やFortranといった、CPUの特性を最大限に引き出す言語で書かれており、現代のほとんどの数値計算ライブラリやアプリケーションの基盤となっています。
Julia は、デフォルトで OpenBLAS などの高性能なBLAS実装を使用しており、ユーザーが意識することなく高速な線形代数演算を享受できます。LinearAlgebra.BLAS
モジュールは、このBLASライブラリと直接やり取りするための機能を提供します。
LinearAlgebra.BLAS
の主な役割と特徴
-
低レベルなBLAS関数のラッパー
LinearAlgebra.BLAS
は、BLASの基本的な関数(例:gemm!
(行列の積),dot
(内積),axpy!
(ベクトルスケーリングと加算) など)をJuliaから呼び出すためのラッパーを提供します。これらの関数は、Float64
,Float32
,Complex128
,Complex64
といった様々な数値型に対応しています。 慣例として、入力配列を上書きする関数は名前に!
が付きます(例:gemm!
)。 -
パフォーマンスの最適化
Juliaが内部的に使用するBLASライブラリ(デフォルトではOpenBLASが多いですが、Intel MKLやApple Accelerateなどに切り替えることも可能です)は、高度に最適化されており、CPUのSIMD命令(単一命令複数データ)やマルチスレッド処理を効果的に利用します。これにより、大規模な行列演算において非常に高いパフォーマンスを実現します。 -
スレッド数の制御
LinearAlgebra.BLAS
を介して、BLASが使用するスレッド数を制御できます。BLAS.get_num_threads()
: 現在のBLASスレッド数を取得します。BLAS.set_num_threads(n)
: BLASが使用するスレッド数をn
に設定します。 これにより、マルチコアCPUでの線形代数演算の並列化を細かく調整できます。
-
使用しているBLASライブラリの確認
BLAS.vendor()
関数を使って、現在JuliaがどのBLASライブラリ(例えば、:openblas
や:mkl
)を使用しているかを確認できます。 -
外部BLASライブラリへの切り替え
Juliaのlibblastrampoline
という仕組みにより、実行時に別のBLASライブラリ(例: MKL.jl や BLISBLAS.jl パッケージを使用してIntel MKLやBLIS)に切り替えることが可能です。これは、特定のハードウェア環境で最適なパフォーマンスを引き出すために役立ちます。
なぜ LinearAlgebra.BLAS
が重要なのか?
- 標準化
BLASは業界標準であるため、異なるプログラミング言語や環境間でのコードの互換性やパフォーマンスの予測可能性を高めます。 - 安定性
BLASルーチンは長年にわたって開発・テストされており、非常に安定しています。 - 速度
大規模な数値計算、特に科学技術計算や機械学習において、行列演算は計算時間の大部分を占めることがよくあります。BLASによる最適化されたルーチンは、これらの計算を劇的に高速化します。
LinearAlgebra.BLAS
は、Juliaの高速な線形代数演算の基盤となる部分ですが、時として問題が発生することもあります。ここでは、よくあるエラーとその解決策について説明します。
BLASスレッド数の設定ミスによるパフォーマンス問題
問題
線形代数演算が期待通りに高速ではない、またはCPUのコアが十分に活用されていないように見える。
原因
BLASが使用するスレッド数が正しく設定されていない可能性があります。特に、シングルスレッドで動作するように設定されている場合や、逆に利用可能なコア数よりもはるかに多いスレッド数が設定されている場合、パフォーマンスが低下することがあります。
トラブルシューティング
- スレッド数の設定
最適なスレッド数は、CPUの物理コア数や計算の規模に依存します。通常は物理コア数と同じか、少し少なめに設定するのが良いとされています。
BLASスレッドは、通常Juliaのスレッドとは独立して設定されます。Juliaのマルチスレッド設定はusing LinearAlgebra BLAS.set_num_threads(Sys.CPU_threads()) # システムの論理コア数に設定 # または、手動で設定 # BLAS.set_num_threads(4) # 例: 4スレッドに設定
JULIA_NUM_THREADS
環境変数で制御されますが、BLASのスレッドはBLAS.set_num_threads()
で制御されるか、使用しているBLASライブラリ(OpenBLAS, MKLなど)の環境変数(例:OPENBLAS_NUM_THREADS
,MKL_NUM_THREADS
)によっても影響を受けます。環境変数とset_num_threads
の両方が設定されている場合、通常は環境変数が優先されます。 - 現在のスレッド数の確認
using LinearAlgebra println("BLAS threads: ", BLAS.get_num_threads())
異なるBLASライブラリのロードによる競合や不安定性
問題
Juliaで線形代数演算を実行すると、クラッシュする、予期せぬ結果になる、または非常に遅くなる。複数のパッケージが異なるBLASライブラリをロードしようとする場合に発生しやすい。
原因
システムに複数のBLASライブラリがインストールされており、Juliaが予期しないBLASライブラリをロードしてしまったり、異なるBLASライブラリ間で競合が発生したりすることがあります。特に、Anaconda環境やMKLを別途インストールしている場合によく見られます。
トラブルシューティング
- 特定のBLASライブラリの使用
Julia 1.7以降では、libblastrampoline
という仕組みにより、Juliaが使用するBLASライブラリをより簡単に制御できます。MKL.jlなどのパッケージを使うことで、Intel MKLを明示的に使用するように設定できます。
ただし、これにより別の問題が発生する場合もあるので、注意が必要です。# MKL.jl をインストール # using Pkg; Pkg.add("MKL") using MKL using LinearAlgebra println("BLAS vendor after MKL: ", BLAS.vendor()) # おそらく :mkl と表示される
- 環境変数の確認
LD_PRELOAD
(Linux) やDYLD_LIBRARY_PATH
(macOS) といった環境変数が、システムにインストールされている別のBLASライブラリをJuliaに強制的にロードさせていないか確認します。これらの環境変数が設定されている場合、一時的に解除してJuliaを起動してみることで問題が解決する場合があります。 - 使用しているBLASライブラリの確認
通常はusing LinearAlgebra println("BLAS vendor: ", BLAS.vendor())
:openblas
と表示されます。:mkl
などが表示される場合は、MKLがロードされていることを意味します。
データ型の不一致や不正な引数
問題
BLAS.gemm!
などの低レベルBLAS関数を直接呼び出した際に、MethodError
やBoundsError
が発生する。
原因
BLAS関数は、特定のデータ型(Float64
, Float32
, Complex128
, Complex64
など)と配列の次元に厳密です。Juliaの型システムがそれを強制するため、間違った型や次元の配列を渡すとエラーになります。また、Fortranの列優先(Column-major)順序とJuliaの列優先順序は一致していますが、引数の意味(例: transpose
フラグ)を誤解している場合もあります。
トラブルシューティング
- 配列の次元とサイズの確認
行列の積の場合、行列の次元が適切に合致しているか確認します。例えば、Am×kと Bk×nの積は Cm×nになります。 - データ型の一致
渡す行列やベクトルの要素型が、BLAS関数がサポートする型と一致していることを確認します。A = rand(Float32, 3, 3) # Float32 の行列 B = rand(Float32, 3, 3) # Float32 の行列 C = zeros(Float32, 3, 3) # 結果格納用も同じ型 LinearAlgebra.BLAS.gemm!('N', 'N', 1.0f0, A, B, 0.0f0, C) # スカラーもFloat32にする
- 関数のシグネチャを確認
?BLAS.gemm!
のようにJuliaのREPLでヘルプを確認し、期待される引数の型と数を確認します。
特定のBLAS実装のバグや非互換性
問題
特定のBLASライブラリ(例: 古いバージョンのOpenBLASやMKL)を使用している場合にのみ、線形代数演算で誤った結果が出たり、クラッシュしたりする。
原因
ごく稀に、特定のBLASライブラリのバージョンにバグが存在したり、Juliaの特定のバージョンとの間で非互換性があったりすることがあります。
トラブルシューティング
- 別のBLASライブラリへの切り替え
可能であれば、OpenBLASからMKLへ、またはその逆へ切り替えてみて、問題が解決するかどうか確認します。これは、MKL.jl
などのパッケージを使用するか、環境変数を設定することで行えます。 - BLASライブラリのバージョン確認
使用しているBLASライブラリのバージョンを確認し、既知のバグがないか、そのライブラリのドキュメントや変更ログを調べます。 - Juliaのアップデート
最新の安定版Juliaにアップデートすることで、BLAS関連の依存関係や内部的な問題が修正されている可能性があります。
環境変数によるBLASの挙動の変化
問題
ある環境では正常に動作するが、別の環境では線形代数演算のパフォーマンスが大きく異なる、またはエラーが発生する。
原因
BLASライブラリの挙動は、特定の環境変数によって制御されることがあります。例えば、OpenBLASはOPENBLAS_NUM_THREADS
、Intel MKLはMKL_NUM_THREADS
などでスレッド数が制御されます。これらの環境変数が異なる環境で異なる値に設定されていると、パフォーマンスや安定性に影響が出ることがあります。
- 明示的な設定
スクリプトの先頭でBLAS.set_num_threads()
を呼び出すことで、環境変数に依存せずにスレッド数を明示的に設定できます。using LinearAlgebra # 環境変数を上書きしたい場合や、確実な設定をしたい場合 BLAS.set_num_threads(Sys.CPU_threads())
- 関連する環境変数の確認
Juliaを起動する前に、OPENBLAS_NUM_THREADS
,MKL_NUM_THREADS
,OMP_NUM_THREADS
などの環境変数がどのように設定されているか確認します。
通常、Juliaで線形代数演算を行う場合、*
演算子や高レベルな関数(例: inv
, det
, eigen
など)を使用することがほとんどです。これらの高レベルな関数は、内部で自動的にLinearAlgebra.BLAS
を介して最適化されたBLASルーチンを呼び出すため、ユーザーが直接BLAS
モジュールの関数を呼び出す必要は稀です。
しかし、LinearAlgebra.BLAS
モジュールの関数を直接呼び出すことで、より低レベルな制御が可能になり、特定の最適化やデバッグ、あるいはBLASの動作を理解するのに役立ちます。
以下に、いくつかの主要なBLAS
関数の例を示します。
BLASスレッド数の設定と確認
BLASが使用するCPUスレッド数を制御できます。これはパフォーマンスチューニングに重要です。
using LinearAlgebra
# 現在のBLASスレッド数を確認
current_threads = BLAS.get_num_threads()
println("現在のBLASスレッド数: ", current_threads)
# BLASスレッド数を変更(例: 4スレッドに設定)
# ご自身のCPUコア数に合わせて調整してください
BLAS.set_num_threads(4)
println("変更後のBLASスレッド数: ", BLAS.get_num_threads())
# システムの論理コア数に設定し直す(推奨される設定の一つ)
BLAS.set_num_threads(Sys.CPU_threads())
println("システムの論理コア数に設定: ", BLAS.get_num_threads())
# 使用しているBLASベンダーの確認
println("使用しているBLASベンダー: ", BLAS.vendor())
出力例
現在のBLASスレッド数: 8
変更後のBLASスレッド数: 4
システムの論理コア数に設定: 8
使用しているBLASベンダー: :openblas
(Sys.CPU_threads()
は論理コア数を返します。私の環境では8コアです。)
ベクトルとベクトルの内積 (BLAS Level 1: dot)
dot
関数は、2つのベクトルの内積を計算します。これはBLAS Level 1の操作です。
using LinearAlgebra
# 浮動小数点ベクトル
x = [1.0, 2.0, 3.0]
y = [4.0, 5.0, 6.0]
# BLASのdot関数を呼び出す
# LinearAlgebra.BLAS.dot(n, x, incx, y, incy)
# n: ベクトルの長さ
# x, y: ベクトル
# incx, incy: ベクトル要素のインクリメント(通常は1)
result_blas = BLAS.dot(length(x), x, 1, y, 1)
println("BLAS.dotによる内積: ", result_blas)
# Juliaの標準dot関数(内部でBLASを使用)
result_julia = dot(x, y)
println("Juliaのdot関数による内積: ", result_julia)
# 結果の確認
@assert result_blas ≈ result_julia
出力
BLAS.dotによる内積: 32.0
Juliaのdot関数による内積: 32.0
ベクトルスケーリングと加算 (BLAS Level 1: axpy!)
axpy!
は、y=αx+y を計算します。!
が付く関数は、引数をインプレースで(元のメモリを上書きして)変更します。
using LinearAlgebra
α = 2.0
x = [1.0, 2.0, 3.0]
y = [10.0, 20.0, 30.0]
y_original = copy(y) # yの元の値を保持
println("元のy: ", y)
# BLASのaxpy!関数を呼び出す
# LinearAlgebra.BLAS.axpy!(n, alpha, x, incx, y, incy)
# n: ベクトルの長さ
# alpha: スカラー定数
# x, y: ベクトル
# incx, incy: インクリメント
BLAS.axpy!(length(x), α, x, 1, y, 1)
println("BLAS.axpy!後のy: ", y)
# 期待される結果(手動計算)
expected_y = α .* x .+ y_original
println("期待されるy: ", expected_y)
@assert y ≈ expected_y
出力
元のy: [10.0, 20.0, 30.0]
BLAS.axpy!後のy: [12.0, 24.0, 36.0]
期待されるy: [12.0, 24.0, 36.0]
行列と行列の積 (BLAS Level 3: gemm!)
gemm!
は、C=αAB+βC を計算する汎用行列積ルーチンです。これはBLAS Level 3の操作であり、非常に重要です。
using LinearAlgebra
# 行列の定義
A = [1.0 2.0; 3.0 4.0] # 2x2 行列
B = [5.0 6.0; 7.0 8.0] # 2x2 行列
C = zeros(2, 2) # 結果を格納する行列 (ゼロで初期化)
# スカラー定数
α = 1.0 # alpha
β = 0.0 # beta (Cを初期化するので0にする)
println("A:\n", A)
println("B:\n", B)
println("初期C:\n", C)
# BLAS.gemm! 関数を呼び出す
# LinearAlgebra.BLAS.gemm!(transA, transB, alpha, A, B, beta, C)
# transA: 'N' (no transpose) または 'T' (transpose) または 'C' (conjugate transpose) for A
# transB: 'N' (no transpose) または 'T' (transpose) または 'C' (conjugate transpose) for B
# alpha: スカラー係数
# A, B: 入力行列
# beta: スカラー係数
# C: 結果を格納する行列(上書きされる)
BLAS.gemm!('N', 'N', α, A, B, β, C)
println("BLAS.gemm!によるA*Bの結果C:\n", C)
# Juliaの標準行列積演算子(内部でBLASを使用)
C_julia = A * B
println("Juliaの*演算子によるA*Bの結果:\n", C_julia)
@assert C ≈ C_julia
出力
A:
[1.0 2.0; 3.0 4.0]
B:
[5.0 6.0; 7.0 8.0]
初期C:
[0.0 0.0; 0.0 0.0]
BLAS.gemm!によるA*Bの結果C:
[19.0 22.0; 43.0 50.0]
Juliaの*演算子によるA*Bの結果:
[19.0 22.0; 43.0 50.0]
gemm!のtransAとtransBについて
'N'
は行列をそのまま使うことを意味します。'T'
は転置して使うことを意味します。
例えば、C=ATB を計算したい場合は、BLAS.gemm!('T', 'N', α, A, B, β, C)
とします。
特定のデータ型での使用例
BLAS関数は、Float64
だけでなく、Float32
やComplexF64
(Complex128
)などの型にも対応しています。
using LinearAlgebra
# Float32 (単精度浮動小数点数) の行列積
A_f32 = Float32[1.0 2.0; 3.0 4.0]
B_f32 = Float32[5.0 6.0; 7.0 8.0]
C_f32 = zeros(Float32, 2, 2)
BLAS.gemm!('N', 'N', 1.0f0, A_f32, B_f32, 0.0f0, C_f32) # スカラーもFloat32にする
println("\nFloat32の行列積の結果C_f32:\n", C_f32)
# ComplexF64 (複素数) の内積
x_cplx = [1.0 + 2.0im, 3.0 + 4.0im]
y_cplx = [5.0 + 6.0im, 7.0 + 8.0im]
# 複素数の内積は共役転置(dotの場合は共役)が含まれる
# BLASのdotcは複素数ベクトルの共役内積 (conjugate dot product) を計算
result_cplx_blas = BLAS.dotc(length(x_cplx), x_cplx, 1, y_cplx, 1)
println("ComplexF64のBLAS.dotcによる内積: ", result_cplx_blas)
# Juliaの標準dot関数も複素数に対応
result_cplx_julia = dot(x_cplx, y_cplx)
println("ComplexF64のJuliaのdot関数による内積: ", result_cplx_julia)
@assert result_cplx_blas ≈ result_cplx_julia
Float32の行列積の結果C_f32:
[19.0f0 22.0f0; 43.0f0 50.0f0]
ComplexF64のBLAS.dotcによる内積: (86.0 - 8.0im)
ComplexF64のJuliaのdot関数による内積: (86.0 - 8.0im)
ここでは、BLASを直接操作する代わりに、Juliaで高速な線形代数演算を行うための一般的な代替方法と、それらの利用シーンについて説明します。
Juliaの標準的な行列演算子 (*)
最も一般的で推奨される方法は、Juliaの標準的な行列乗算演算子である *
を使用することです。これは、内部でLinearAlgebra.BLAS.gemm!
などのBLAS Level 3ルーチンを自動的に呼び出すため、最高のパフォーマンスを提供します。
例
using LinearAlgebra
A = rand(1000, 1000) # 1000x1000のランダム行列
B = rand(1000, 1000)
# 最も一般的で効率的な行列積の計算方法
C = A * B
# ベクトルと行列の積も同様
v = rand(1000)
w = A * v
利点
- 型推論の恩恵
Juliaのコンパイラが型情報を最大限に活用し、効率的なコードを生成します。 - 最高のパフォーマンス
内部で最適化されたBLASルーチンを自動的に利用します。ユーザーが意識する必要がありません。 - シンプルで読みやすいコード
数学的な表記に近く、直感的です。
いつ使うか
線形代数演算を行うほとんどすべてのケースで、この方法を使用すべきです。
LinearAlgebra モジュールの高レベル関数
LinearAlgebra
モジュールは、行列式(det
), 逆行列(inv
), 固有値分解(eigen
), 特異値分解(svd
)など、様々な線形代数操作のための高レベルな関数を提供します。これらの関数も、内部でBLASやLAPACK(Linear Algebra PACKage, BLASの上に構築されたより複雑な線形代数計算ライブラリ)の最適化されたルーチンを呼び出します。
例
using LinearAlgebra
M = rand(5, 5)
# 行列式
d = det(M)
println("行列式: ", d)
# 逆行列
Minv = inv(M)
println("逆行列:\n", Minv)
# 固有値分解
# E.values は固有値のベクトル、E.vectors は固有ベクトルを行列にしたもの
E = eigen(M)
println("固有値: ", E.values)
println("固有ベクトル:\n", E.vectors)
# ベクトル内積(dotも内部でBLAS Level 1を使用)
x = rand(10)
y = rand(10)
dot_product = dot(x, y)
println("内積: ", dot_product)
利点
- 高性能
内部でBLAS/LAPACKを自動的に利用します。 - 信頼性
長年にわたってテストされた、堅牢なアルゴリズムに基づいています。 - 機能が豊富
線形代数の幅広い問題に対応できます。
いつ使うか
特定の線形代数問題を解く必要があるが、低レベルな実装の詳細に立ち入る必要がない場合。
StaticArrays.jl (小さな配列/行列向け)
非常に小さな固定サイズの配列や行列(例: 2x2, 3x3, 4x4など)を扱う場合、StaticArrays.jl
パッケージは、BLASを使うよりもはるかに高速な演算を提供できます。これは、コンパイル時に配列のサイズが既知であるため、メモリ割り当てを避け、ループ展開などのコンパイル時最適化を最大限に活用できるためです。
例
using StaticArrays, BenchmarkTools
# 静的配列の定義 (コンパイル時にサイズが決定)
A_s = @SMatrix [1.0 2.0; 3.0 4.0]
B_s = @SMatrix [5.0 6.0; 7.0 8.0]
# 静的配列の積は、BLASを介さずコンパイル時最適化で実行される
C_s = A_s * B_s
println("StaticArraysによる行列積:\n", C_s)
# パフォーマンス比較(小規模な行列の場合)
# 動的配列 (通常のArray)
A_d = rand(2, 2)
B_d = rand(2, 2)
@btime $A_d * $B_d # $ を付けてグローバル変数からの脱出を防ぐ
# 静的配列
A_s = @SMatrix rand(2, 2)
B_s = @SMatrix rand(2, 2)
@btime $A_s * $B_s
利点
- ゼロアロケーション
メモリ割り当てが発生しないため、ガベージコレクションのオーバーヘッドがありません。 - 極めて高いパフォーマンス
特に小さなサイズの場合、BLASを上回る速度が出ることがあります。
いつ使うか
シミュレーションなどで、多くの小さな線形代数演算を繰り返す場合(例: 物理エンジンの剛体変換、グラフィックスの行列変換)。
LoopVectorization.jl (手動最適化)
LoopVectorization.jl
は、手書きのループで線形代数演算を実装する際に、SIMD(Single Instruction, Multiple Data)命令やマルチスレッド化を自動的に適用し、大幅な高速化を図るためのパッケージです。これは、BLASが提供しない特定の操作や、BLASのオーバーヘッドを避けたい場合に有用です。
例
(簡単な行列積のループを手動でベクトル化)
using LoopVectorization
using BenchmarkTools
function my_matrix_multiply!(C, A, B)
@turbo for i in axes(A, 1), j in axes(B, 2)
Cij = zero(eltype(C))
for k in axes(A, 2)
Cij += A[i, k] * B[k, j]
end
C[i, j] = Cij
end
return C
end
A = rand(100, 100)
B = rand(100, 100)
C_custom = zeros(100, 100)
C_blas = zeros(100, 100)
@btime my_matrix_multiply!($C_custom, $A, $B) # 手動ループ + @turbo
@btime mul!($C_blas, $A, $B) # BLASによる行列積(推奨)
利点
- 強力な最適化
コンパイラによるベクトル化や並列化を強化します。 - 柔軟性
BLASが提供しないカスタムな線形代数操作を最適化できます。
いつ使うか
- 学習目的や、Juliaのパフォーマンス最適化の仕組みを理解するため。
- 既存のBLASルーチンよりもさらに細かい制御やチューニングが必要な場合(ただし、これは稀です)。
- BLASが提供しない特殊な行列操作を実装する場合。
特定のBLAS/LAPACK実装 (MKL.jl, BLISBLAS.jlなど)
JuliaのLinearAlgebra
はデフォルトでOpenBLASを使用しますが、パフォーマンスが重要な場合は、より高速なBLAS実装(例えばIntel MKLやBLIS)に切り替えることを検討できます。これはMKL.jl
やBLISBLAS.jl
といったパッケージを通じて行われます。
例
# まず MKL.jl をインストールします
# using Pkg
# Pkg.add("MKL")
using MKL
using LinearAlgebra
# MKLがロードされたことを確認
println("現在のBLASベンダー: ", BLAS.vendor())
A = rand(1000, 1000)
B = rand(1000, 1000)
# MKLが有効になっている場合、この計算はMKLのgemmによって実行される
C = A * B
利点
- 安定性
企業レベルで開発・テストされたライブラリです。 - 最高レベルのパフォーマンス
特定のハードウェア(Intel CPU上のMKLなど)で、最高のパフォーマンスを発揮することがあります。
いつ使うか
- 特定のハードウェア(Intel CPUなど)で最大限のパフォーマンスを引き出したい場合。
- 大規模な数値シミュレーションや機械学習モデルのトレーニングなど、線形代数計算がボトルネックになっている場合。
Juliaで線形代数プログラミングを行う際の「代替方法」は、実際にはJuliaの組み込みの強力な機能や、高性能なパッケージを活用することに他なりません。
- 究極のパフォーマンス
MKL.jl
などのパッケージでBLAS実装を切り替えることを検討します。 - 手動最適化/カスタム操作
LoopVectorization.jl
は、より低レベルなループ最適化を提供します。 - 非常に小さい配列/行列
StaticArrays.jl
は、BLASよりも高速になる可能性があります。 - ほとんどの場合
Juliaの標準的な演算子 (*
) やLinearAlgebra
モジュールの高レベル関数を使用するのが最善です。これらは内部でBLASを効率的に利用します。