Juliaの線形代数:BLAS.gemv()のエラー対処とトラブルシューティング

2025-05-26

BLAS とは何か?

まず、BLAS (Basic Linear Algebra Subprograms) について簡単に説明します。BLAS は、線形代数における基本的な操作(ベクトルとベクトルの積、行列とベクトルの積、行列と行列の積など)を効率的に行うための標準的なルーチンの集合です。多くの数値計算ライブラリや科学技術計算ソフトウェアは、内部で BLAS を利用することで、高いパフォーマンスを実現しています。

BLAS のルーチンは、その計算量に基づいて「レベル」に分類されています。

  • Level 3 BLAS: 行列-行列演算 (O(n3))
  • Level 2 BLAS: 行列-ベクトル演算 (O(n2))
  • Level 1 BLAS: ベクトル-ベクトル演算 (O(n))

gemvLevel 2 BLAS に分類されます。

gemv とは何か?

gemv は "general matrix-vector multiplication" の略で、以下の形式の行列-ベクトル積を計算するルーチンです。

y=αAx+βy

ここで、

  • β: スカラー係数
  • α: スカラー係数
  • y: ベクトル
  • x: ベクトル
  • A: 一般行列 (general matrix)

です。

Juliaにおける LinearAlgebra.BLAS.gemv()

Juliaでは、標準ライブラリの LinearAlgebra モジュールを通じて BLAS ルーチンにアクセスできます。LinearAlgebra.BLAS.gemv() は、この BLASgemv ルーチンを直接呼び出すための関数です。

LinearAlgebra.BLAS.gemv! の基本的な使い方

Juliaでは、多くの線形代数関数に破壊的 (in-place) なバージョンがあり、関数名の末尾に ! が付きます。gemv! もその一つで、結果を既存のベクトル y に上書きします。

最も一般的な形式は以下のようになります(実際には引数がもっと多いですが、理解のために簡略化しています)。

LinearAlgebra.BLAS.gemv!(trans, alpha, A, x, beta, y)

各引数の意味は以下の通りです。

  • y: ベクトル y。計算結果がこの y に格納されます。
  • beta: スカラー係数 β。
  • x: ベクトル x。
  • A: 行列。
  • alpha: スカラー係数 α。
  • trans: 行列 A を転置するかどうかを指定します。
    • 'N' (No transpose): A をそのまま使用します。
    • 'T' (Transpose): A を転置して使用します (AT)。
    • 'C' (Conjugate transpose): A を共役転置して使用します (AH)。複素数の場合に用いられます。

具体例

簡単な例を見てみましょう。

using LinearAlgebra

A = [1.0 2.0; 3.0 4.0]
x = [1.0; 2.0]
y = [0.0; 0.0] # 結果を格納するベクトルを初期化

# y = 1.0 * A * x + 0.0 * y  つまり y = A * x を計算
LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y)
println("y = A * x: ", y) # 出力: y = A * x: [5.0, 11.0]

# y = 2.0 * A * x + 1.0 * y  つまり y = 2 * (A * x) + y_old を計算
# (y_old は上の計算結果 [5.0, 11.0])
A = [1.0 2.0; 3.0 4.0]
x = [1.0; 2.0]
y = [5.0; 11.0] # 前の結果を使用

LinearAlgebra.BLAS.gemv!('N', 2.0, A, x, 1.0, y)
println("y = 2 * A * x + y_old: ", y) # 出力: y = 2 * A * x + y_old: [15.0, 33.0]

通常、Juliaで行列-ベクトル積を計算する場合、A * x のようにシンプルな記法を使用します。Juliaの内部では、多くの場合、この記法が最適な BLAS ルーチン(この場合は gemv)に自動的にディスパッチされます。

しかし、以下のような場合に LinearAlgebra.BLAS.gemv() を直接使うことが考えられます。

  1. パフォーマンスの最適化: 非常に低レベルな部分で、特定の BLAS の挙動を細かく制御したい場合。例えば、alphabeta を使って中間結果を蓄積するような複雑な計算グラフを構築する際など。
  2. 既存のC/Fortranコードとの連携: BLAS はCやFortranで実装されており、既存のコードベースをJuliaに移植する際に、BLAS インターフェースを直接呼び出す必要が生じることがあります。
  3. BLASの実装固有の機能へのアクセス: ごく稀に、特定の BLAS 実装(OpenBLAS, MKLなど)が提供する、標準 BLAS 仕様にはない、しかし gemv 関連の低レベルな機能にアクセスしたい場合。

注意点: ほとんどのユーザーは、通常 A * x のようにJuliaの直感的な構文を使うべきです。Juliaの線形代数ライブラリは、内部で最適な BLAS 関数を呼び出すように設計されており、ほとんどの場合、直接 BLAS.gemv! を呼び出すよりも効率的または同等のパフォーマンスが得られます。LinearAlgebra.gemv! (非破壊的バージョン) や mul! 関数(汎用的な行列積関数)も、内部で BLAS.gemv! を利用することがあります。



gemv! は、y = α A x + β y の計算を行います。この関数を呼び出す際に最もよくある問題は、引数の型、次元、または順序の不一致です。

引数の型が間違っている (MethodError)

BLAS.gemv! は、特定の数値型(通常はFloat32, Float64, ComplexF32, ComplexF64)に対して最適化されています。異なる型の引数を与えるとMethodErrorが発生する可能性があります。

エラー例

using LinearAlgebra

A = [1 2; 3 4] # Int型の行列
x = [1; 2]     # Int型のベクトル
y = [0; 0]

LinearAlgebra.BLAS.gemv!('N', 1, A, x, 0, y)

考えられるエラーメッセージ
MethodError: no method matching gemv!(::Char, ::Int64, ::Matrix{Int64}, ::Vector{Int64}, ::Int64, ::Vector{Int64})

原因
A, x, y の要素が Int 型であるため、gemv! が対応するメソッドを見つけられません。BLASルーチンは浮動小数点数(または複素数)の計算に特化しています。

トラブルシューティング
引数を適切な浮動小数点型に変換します。

using LinearAlgebra

A = [1.0 2.0; 3.0 4.0] # Float64型の行列
x = [1.0; 2.0]         # Float64型のベクトル
y = [0.0; 0.0]

LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y)
println(y)

または、. をつけてリテラルを浮動小数点数にします。

A = [1.0 2.0; 3.0 4.0]
x = [1.0; 2.0]
y = zeros(2) # zerosでFloat64のベクトルを生成
LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y)

引数の次元が合っていない (DimensionMismatch)

行列-ベクトル積の計算では、行列とベクトルの次元が一致している必要があります。

  • Aがm×n 行列の場合:
    • trans = 'N' の場合 (A * x):x は n 要素のベクトル、y は m 要素のベクトルである必要があります。
    • trans = 'T' または 'C' の場合 (Aᵀ * x または Aᴴ * x):x は m 要素のベクトル、y は n 要素のベクトルである必要があります。

エラー例

using LinearAlgebra

A = rand(3, 2) # 3x2行列
x = rand(3)    # 3要素ベクトル (nが2なので合わない)
y = zeros(3)

LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y)

考えられるエラーメッセージ
DimensionMismatch("matrix A has dimensions (3, 2) but vector x has length 3") のように、具体的な次元の不一致を指摘するメッセージが表示されます。

原因
上記の例では、A は 3×2 行列ですが、x は長さ3のベクトルです。A * x を計算する場合、xA の列数(この場合2)と同じ長さである必要があります。

トラブルシューティング
A の次元に合わせて xy の次元を調整します。

using LinearAlgebra

A = rand(3, 2) # 3x2行列
x = rand(2)    # 長さ2のベクトル
y = zeros(3)   # 長さ3のベクトル

LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y)
println(y)

# 転置の場合の例 ('T')
A = rand(3, 2) # 3x2行列
x_T = rand(3)  # Aᵀの行数(Aの列数)と同じ長さ
y_T = zeros(2) # Aᵀの列数(Aの行数)と同じ長さ

LinearAlgebra.BLAS.gemv!('T', 1.0, A, x_T, 0.0, y_T)
println(y_T)

引数の順序や数が間違っている (MethodError)

gemv! の引数の順序は固定されています。間違った順序で引数を渡したり、引数の数が不足していると MethodError になります。

エラー例

using LinearAlgebra

A = rand(2,2)
x = rand(2)
y = zeros(2)

# alphaとAの順序を間違える
LinearAlgebra.BLAS.gemv!('N', A, 1.0, x, 0.0, y)

考えられるエラーメッセージ
これも MethodError となりますが、引数の型と位置から、どのメソッドにもマッチしないことが示されます。

原因
関数のシグネチャと異なる引数の渡し方をしている。

トラブルシューティング
JuliaのドキュメントやREPLのヘルプ(?LinearAlgebra.BLAS.gemv!)で正しい引数順序を確認し、それに従って引数を渡します。

trans 引数の値が不正 (ArgumentError)

trans 引数は 'N', 'T', 'C' のいずれかである必要があります。これ以外の文字を指定するとエラーになります。

エラー例

using LinearAlgebra

A = rand(2,2)
x = rand(2)
y = zeros(2)

LinearAlgebra.BLAS.gemv!('X', 1.0, A, x, 0.0, y)

考えられるエラーメッセージ
ArgumentError: trans must be 'N', 'T', or 'C' のようなメッセージが表示されます。

原因
無効な転置オプションが指定されている。

トラブルシューティング
trans 引数を正しい値に修正します。

配列のメモリーレイアウト(Juliaではあまり問題にならないがBLASの概念として)

BLASルーチンは通常、Fortranの列優先(Column-Major)配列レイアウトを想定しています。Juliaの配列はデフォルトで列優先なので、通常はこの点で問題は発生しません。しかし、もしCライブラリとのFFI(Foreign Function Interface)などで、行優先の配列を渡そうとすると、予期しない結果になる可能性があります。

トラブルシューティング
Juliaの配列はデフォルトで列優先なので、通常は意識する必要がありません。もしCなど他の言語で作られた行優先の配列データをJuliaに渡す場合は、転置を考慮したり、データコピーを行うなどの工夫が必要です。しかし、LinearAlgebra.BLAS.gemv!()をJulia内で直接使う限り、この問題に遭遇することは稀です。

  • ?LinearAlgebra.BLAS.gemv! をREPLで実行すると、関数の完全なシグネチャとドキュメンテーションが表示され、正しい引数を確認できます。
  • トラブルシューティングの第一歩は、エラーメッセージを注意深く読むことです。 Juliaのエラーメッセージは非常に具体的で、問題の原因を特定するのに役立ちます。
  • ほとんどの場合、A * x (非破壊的) または mul!(y, A, x) (破壊的) のように、より高レベルで抽象化されたJuliaの構文を使用すべきです。これらは内部で最適なBLASルーチンを呼び出し、引数の型チェックや次元チェックも適切に行ってくれます。
  • LinearAlgebra.BLAS.gemv!() は、Juliaの線形代数計算の低レベルな詳細であり、通常は直接使用する必要はありません


繰り返しになりますが、Julia では通常 A * xmul!(y, A, x) のような高レベルな関数を使用し、BLAS.gemv!() を直接呼び出すことは稀です。しかし、この低レベルな関数がどのように動作するかを理解することは、Julia の線形代数処理の内部を知る上で役立ちます。

例1: 最も基本的な gemv! の使用方法 (y=Ax)

最も単純な行列-ベクトル積 $y = Ax$ を計算します。この場合、alpha1.0beta0.0 と設定します。trans'N' (No transpose) とします。

using LinearAlgebra

# (1) 行列 A とベクトル x を定義
A = [1.0 2.0;
     3.0 4.0;
     5.0 6.0] # 3x2 行列

x = [1.0;
     2.0]    # 2要素ベクトル

# (2) 結果を格納するためのベクトル y を初期化
# Aは3x2、xは2x1なので、A*xの結果は3x1になる。
# yの要素の型はA,xに合わせてFloat64にする。
y = zeros(3) # [0.0, 0.0, 0.0]

println("--- 例1: y = A * x の計算 ---")
println("A:\n", A)
println("x: ", x)
println("初期の y: ", y)

# (3) BLAS.gemv! を呼び出し
# gemv!(trans, alpha, A, x, beta, y)
LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y)

println("計算後の y: ", y)

# (4) 比較のために通常の行列積も計算
y_standard = A * x
println("A * x (標準): ", y_standard)

# 結果は一致することを確認
@assert y ≈ y_standard
println("結果は標準の計算と一致します。\n")

# 計算の内訳:
# y[1] = (1.0 * 1.0) + (2.0 * 2.0) = 1.0 + 4.0 = 5.0
# y[2] = (3.0 * 1.0) + (4.0 * 2.0) = 3.0 + 8.0 = 11.0
# y[3] = (5.0 * 1.0) + (6.0 * 2.0) = 5.0 + 12.0 = 17.0

例2: 転置行列との積 (y=ATx)

trans 引数を 'T' に設定することで、行列 A の転置 AT とベクトル x の積を計算できます。

using LinearAlgebra

# (1) 行列 A とベクトル x を定義
A = [1.0 2.0 3.0;
     4.0 5.0 6.0] # 2x3 行列

x = [1.0;
     2.0]        # 2要素ベクトル (A^Tの行数と同じ)

# (2) 結果を格納するためのベクトル y を初期化
# A^Tは3x2、xは2x1なので、A^T*xの結果は3x1になる。
y = zeros(3)

println("--- 例2: y = A^T * x の計算 ---")
println("A:\n", A)
println("x: ", x)
println("初期の y: ", y)

# (3) BLAS.gemv! を呼び出し ('T' で転置を指定)
LinearAlgebra.BLAS.gemv!('T', 1.0, A, x, 0.0, y)

println("計算後の y: ", y)

# (4) 比較のために通常の転置行列積も計算
y_standard = A' * x
println("A' * x (標準): ", y_standard)

@assert y ≈ y_standard
println("結果は標準の計算と一致します。\n")

# 計算の内訳:
# A^T = [1.0 4.0;
#        2.0 5.0;
#        3.0 6.0]
# y[1] = (1.0 * 1.0) + (4.0 * 2.0) = 1.0 + 8.0 = 9.0
# y[2] = (2.0 * 1.0) + (5.0 * 2.0) = 2.0 + 10.0 = 12.0
# y[3] = (3.0 * 1.0) + (6.0 * 2.0) = 3.0 + 12.0 = 15.0

例3: スカラー係数 alphabeta の使用 (y=αAx+βy)

alphabeta を非ゼロに設定することで、より一般的な線形変換を計算できます。

using LinearAlgebra

# (1) 行列 A とベクトル x を定義
A = [1.0 2.0;
     3.0 4.0] # 2x2 行列

x = [1.0;
     2.0]    # 2要素ベクトル

# (2) 結果を格納するためのベクトル y を初期化 (既存の値を持つ)
y = [10.0;
     20.0]

alpha = 2.0
beta = 0.5

println("--- 例3: y = α A x + β y の計算 ---")
println("A:\n", A)
println("x: ", x)
println("初期の y: ", y)
println("alpha: ", alpha)
println("beta: ", beta)

# (3) BLAS.gemv! を呼び出し
# y = 2.0 * (A * x) + 0.5 * y_old
LinearAlgebra.BLAS.gemv!('N', alpha, A, x, beta, y)

println("計算後の y: ", y)

# (4) 比較のために手動で計算
# まず A*x を計算
Ax = A * x # [5.0, 11.0]

# 次に alpha * Ax を計算
alpha_Ax = alpha * Ax # [10.0, 22.0]

# 最後に beta * y_old を計算 (y_old は [10.0, 20.0])
beta_y_old = beta * [10.0, 20.0] # [5.0, 10.0]

# 最終結果は alpha_Ax + beta_y_old
y_expected = alpha_Ax + beta_y_old # [10.0 + 5.0, 22.0 + 10.0] = [15.0, 32.0]
println("手動計算での y: ", y_expected)

@assert y ≈ y_expected
println("結果は手動計算と一致します。\n")

複素数の場合、trans'C' (共役転置) を指定できます。

using LinearAlgebra

# (1) 複素数行列 A とベクトル x を定義
A = [1.0 + 2.0im  3.0 + 4.0im;
     5.0 + 6.0im  7.0 + 8.0im] # 2x2 複素数行列

x = [1.0 + 1.0im;
     2.0 + 2.0im]              # 2要素複素数ベクトル

# (2) 結果を格納するための複素数ベクトル y を初期化
y = zeros(ComplexF64, 2) # ComplexF64型のベクトル

println("--- 例4: 複素数に対する y = A^H * x の計算 ---")
println("A:\n", A)
println("x: ", x)
println("初期の y: ", y)

# (3) BLAS.gemv! を呼び出し ('C' で共役転置を指定)
# alpha = 1.0, beta = 0.0 で y = A^H * x
LinearAlgebra.BLAS.gemv!('C', 1.0 + 0.0im, A, x, 0.0 + 0.0im, y)

println("計算後の y: ", y)

# (4) 比較のために通常の共役転置行列積も計算
y_standard = A' * x # A' は共役転置を意味する
println("A' * x (標準): ", y_standard)

@assert y ≈ y_standard
println("結果は標準の計算と一致します。\n")

# 計算の内訳:
# A^H = [1.0 - 2.0im  5.0 - 6.0im;
#        3.0 - 4.0im  7.0 - 8.0im]
#
# y[1] = (1.0 - 2.0im) * (1.0 + 1.0im) + (5.0 - 6.0im) * (2.0 + 2.0im)
#      = (1.0 + 1.0im - 2.0im + 2.0) + (10.0 + 10.0im - 12.0im + 12.0)
#      = (3.0 - 1.0im) + (22.0 - 2.0im) = 25.0 - 3.0im
#
# y[2] = (3.0 - 4.0im) * (1.0 + 1.0im) + (7.0 - 8.0im) * (2.0 + 2.0im)
#      = (3.0 + 3.0im - 4.0im + 4.0) + (14.0 + 14.0im - 16.0im + 16.0)
#      = (7.0 - 1.0im) + (30.0 - 2.0im) = 37.0 - 3.0im

これらの例は、LinearAlgebra.BLAS.gemv!()alpha, A, x, beta, ytrans オプションを使って、線形代数の基本的な操作 y = α A x + β y をいかに直接的に実行するかを示しています。

実際のJuliaプログラミングでは、ほとんどの場合、以下のような高レベルな関数や記法を使用します。

  • 破壊的(結果を既存の変数に上書き):
    • mul!(y, A, x) (これは alpha=1.0, beta=0.0gemv! に相当)
    • mul!(y, A, x, alpha, beta) (これは gemv! と完全に同じ機能)
  • 非破壊的(結果を新しい変数に格納):
    • y = A * x
    • y = A' * x (共役転置の場合、実数の場合は転置)
    • y = alpha * A * x + beta * y_old


LinearAlgebra.BLAS.gemv() の代替方法

gemvの主要な機能は行列とベクトルの積を計算することです。Juliaの一般的な記法や関数はこれをより使いやすく、破壊的(in-place)か非破壊的(out-of-place)かを選択できるようにしています。

非破壊的な行列-ベクトル積 (Out-of-place)

これは最も一般的で直感的な方法です。新しいベクトルを生成して結果を格納します。

  • y = alpha * A * x + beta * y_old のような複合演算: Juliaの記法は非常に柔軟なので、alphabeta を含む一般的な線形変換も直接記述できます。

    using LinearAlgebra
    
    A = [1.0 2.0; 3.0 4.0]
    x = [1.0; 2.0]
    y_old = [10.0; 20.0]
    
    alpha = 2.0
    beta = 0.5
    
    # y = alpha * A * x + beta * y_old を計算
    # これは gemv!('N', alpha, A, x, beta, y) と y が新しい変数である点で異なる
    # gemv!はyを上書きするが、この記法は新しいyを生成する
    y_new = alpha * A * x + beta * y_old
    println("y_new = alpha * A * x + beta * y_old: ", y_new)
    
  • y = A' * x: 行列 A の転置(実数の場合)または共役転置(複素数の場合)とベクトル x の積を計算します。これも内部で適切なBLASルーチン(gemv!('T', ...) または gemv!('C', ...) に相当)が呼び出されます。

    using LinearAlgebra
    
    A = [1.0 2.0 3.0; 4.0 5.0 6.0] # 2x3 行列
    x = [1.0; 2.0]                 # 2要素ベクトル
    
    # y = A' * x を計算 (gemv!('T', 1.0, A, x, 0.0, y) に相当)
    y_transposed = A' * x
    println("y = A' * x: ", y_transposed)
    
    # 複素数の例
    A_complex = [1.0+im 2.0-im; 3.0im 4.0+2im]
    x_complex = [1.0-im; 2.0+im]
    
    # y = A^H * x を計算 (gemv!('C', 1.0, A_complex, x_complex, 0.0, y) に相当)
    y_conjugate_transposed = A_complex' * x_complex
    println("y = A^H * x (complex): ", y_conjugate_transposed)
    
  • y = A * x: 最も直接的な行列-ベクトル積の計算です。Juliaのコンパイラは、これが最適なBLASのgemvルーチン(alpha=1.0, beta=0.0 のケース)にディスパッチされるように最適化されています。

    using LinearAlgebra
    
    A = [1.0 2.0; 3.0 4.0]
    x = [1.0; 2.0]
    
    # y = A * x を計算 (gemv!('N', 1.0, A, x, 0.0, y) に相当)
    y = A * x
    println("y = A * x: ", y)
    

破壊的な行列-ベクトル積 (In-place)

結果を既存のベクトルに上書きしたい場合(メモリ割り当てを避けたい場合など)は、LinearAlgebra モジュールの mul! 関数を使用します。

  • mul!(y, A, x, alpha, beta): これが LinearAlgebra.BLAS.gemv!(trans, alpha, A, x, beta, y) の最も直接的な代替です。trans 引数がないことに注意してください。これは行列 A 自体を転置して渡すことで対応します(例: mul!(y, A', x, alpha, beta))。

    using LinearAlgebra
    
    A = [1.0 2.0; 3.0 4.0]
    x = [1.0; 2.0]
    y = [10.0; 20.0] # 初期値を持つy
    
    alpha = 2.0
    beta = 0.5
    
    println("--- 破壊的な計算: mul!(y, A, x, alpha, beta) ---")
    println("初期の y: ", y)
    println("alpha: ", alpha)
    println("beta: ", beta)
    
    # y = alpha * (A * x) + beta * y を計算し、yに上書き
    mul!(y, A, x, alpha, beta)
    println("計算後の y: ", y)
    
    # 転置の例: y = alpha * (A' * x) + beta * y
    A_T = [1.0 2.0 3.0; 4.0 5.0 6.0] # 2x3 行列
    x_T = [1.0; 2.0]                 # 2要素ベクトル
    y_T = [100.0; 200.0; 300.0]      # 3要素ベクトル
    
    println("\n--- 破壊的な計算 (転置): mul!(y, A', x, alpha, beta) ---")
    println("A_T:\n", A_T)
    println("x_T: ", x_T)
    println("初期の y_T: ", y_T)
    
    mul!(y_T, A_T', x_T, alpha, beta) # A'を渡すことで転置を処理
    println("計算後の y_T: ", y_T)
    
  • mul!(y, A, x): ベクトル yA * x の結果を上書きします。これは LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, y) と同等です。

    using LinearAlgebra
    
    A = [1.0 2.0; 3.0 4.0]
    x = [1.0; 2.0]
    y = zeros(2) # 結果を格納する既存のベクトル
    
    println("--- 破壊的な計算: mul!(y, A, x) ---")
    println("初期の y: ", y)
    
    mul!(y, A, x) # y を A * x で上書き
    println("計算後の y: ", y)
    
  1. 可読性と直感性: A * xmul!(y, A, x) は数学的な記法や一般的なプログラミングの慣習に近く、コードの意図が明確です。
  2. 安全性: 高レベルな関数は、引数の型や次元の不一致などを自動的にチェックし、分かりやすいエラーメッセージを提供します。BLAS.gemv!() を直接使う場合、これらのチェックはBLASルーチン自体に依存するため、より低レベルなエラーになりがちです。
  3. 最適化: Juliaのコンパイラと標準ライブラリは、これらの高レベルな関数を内部で最適なBLASルーチンに自動的にディスパッチするように設計されています。通常、直接 BLAS.gemv!() を呼び出すよりも、パフォーマンス上のメリットがあるか、少なくとも同等のパフォーマンスが得られます。JuliaはBLASの実装(OpenBLAS, MKLなど)も自動的に選択・利用します。
  4. 汎用性: mul! は行列-ベクトル積だけでなく、行列-行列積 (mul!(C, A, B)) にも使用できるなど、より広範な線形代数演算に適用できる汎用的なインターフェースを提供します。
  5. 型の抽象化: A * xmul! は、AbstractMatrixAbstractVector などの抽象型で動作するため、具体的な行列やベクトルの種類(Matrix, SparseMatrixCSC, Adjoint, SubArray など)によらず、同じコードで処理できます。これは BLAS.gemv!() が特定の具象型に特化しているのと対照的です。