数値計算の精度向上: Julia `trrfs!()`と誤差解析

2025-05-27

LinearAlgebra.LAPACK.trrfs!()とは何か

JuliaのLinearAlgebraモジュールは、線形代数に関する様々な機能を提供しています。その中でもLAPACKサブモジュールは、高性能な線形代数計算ライブラリであるLAPACK(Linear Algebra PACKage)の低レベルな関数を直接呼び出すためのインターフェースを提供しています。

trrfsは、LAPACKにおける三角形行列の連立一次方程式の解の改良反復 (iterative refinement)誤差推定 (error estimation) を行うルーチンの名前です。

つまり、LinearAlgebra.LAPACK.trrfs!(...)は、以下の目的で使用されます。

  1. 三角形行列の連立一次方程式の解の精度を高める: 既に計算された解に対して、さらに精度を上げるための反復改良を行います。
  2. 解の誤差(前方誤差および後方残差)を推定する: 計算された解がどれくらいの精度を持っているか、その誤差の上限を評価します。

関数名の末尾に!が付いているのは、Juliaの慣習で、この関数が引数として渡された配列の内容を破壊的に変更する(in-place operation) ことを示しています。

具体的な機能と使用場面

trrfsは、以下のような状況で役立ちます。

  • 既知の解の検証: 特定のアルゴリズムで得られた解が、本当に正しい解に近いのかを検証する際に、その誤差を推定するために利用できます。
  • 解の信頼性保証: 金融計算や科学技術計算など、解の精度が非常に重要となる分野で、計算結果が許容できる誤差範囲内にあることを確認するために使用されます。
  • 数値的に不安定な問題: 特定の条件(例えば、行列がほぼ特異である場合など)では、連立一次方程式の解は非常に不安定になり、浮動小数点計算の誤差が大きく影響することがあります。trrfsを使うことで、その解の信頼性を評価したり、より正確な解を得るための改良を試みたりできます。

LAPACKのtrrfsルーチンは、通常、以下のような情報を入力として受け取り、出力として提供します。

入力

  • IWORK: 整数作業用配列。
  • WORK: 作業用配列。
  • BERR: 解の後方残差 (backward error) の推定値(出力用配列)。
  • FERR: 解の前方誤差 (forward error) の推定値(出力用配列)。
  • X: 既に計算されている解 X 。
  • B: 右辺ベクトル B 。
  • A: 三角行列 A 。
  • N: 行列の次数。
  • DIAG: 行列の対角成分が単位行列('U')か非単位行列('N')かを示す文字。
  • TRANS: 連立一次方程式が AX=B、ATX=B、AHX=B のどれであるかを示す文字。
  • UPLO: 三角行列が上三角('U')か下三角('L')かを示す文字。

出力

  • INFO: 処理の成功/失敗を示す整数コード。
  • BERR: 各右辺ベクトルに対応する解の後方残差の推定値。
  • FERR: 各右辺ベクトルに対応する解の前方誤差の推定値。

JuliaのLinearAlgebra.LAPACK.trrfs!()も、これらのLAPACKルーチンの引数をJuliaの型(Matrix{Float64}など)にマッピングして呼び出します。具体的な引数の順序や型については、Juliaの公式ドキュメントやLAPACKのドキュメントを参照するのが最も確実です。



LinearAlgebra.LAPACK.trrfs!()の一般的なエラーとトラブルシューティング

LAPACKException(info)

これはLAPACK関数がエラーを返した場合に最も頻繁に発生するエラーです。infoは整数値で、エラーの種類を示します。

  • info > 0: 計算が収束しなかった、または行列が特異であるなど、数値的な問題が発生したことを意味します。trrfsの場合、これは通常、解の反復改良が所定の回数で収束しなかったことを示します。
    • トラブルシューティング:
      • 入力行列 A が特異またはほぼ特異である: 三角行列 A が特異(非可逆)である場合、厳密な解は存在せず、LAPACKルーチンはエラーを返すか、不安定な結果を生成します。数値計算においては、「ほぼ特異」な行列も同様の問題を引き起こします。
        • A の対角要素に非常に小さい値やゼロがないか確認してください。
        • 可能であれば、問題の条件数(condition number)を調べてください。非常に大きい条件数は、問題が数値的に不安定であることを示唆します。cond(A)で確認できます。
      • 解の精度要求が高すぎる: trrfsは解の精度を改善しようとしますが、元々の問題が非常に悪条件である場合、いくら反復しても望む精度に到達できないことがあります。
      • 数値的な不安定性: 浮動小数点数の丸め誤差が蓄積し、計算が不安定になることがあります。特に、大規模なシステムや非常に小さな値と大きな値が混在するシステムで発生しやすいです。
        • データ型をFloat32からFloat64に、またはComplexF32からComplexF64に変更することで、精度が改善される場合があります。
        • 問題のスケールを調整することを検討してください。
      • FERRBERRの異常値: trrfsの出力である前方誤差(FERR)や後方残差(BERR)がNaN(Not a Number)やInf(Infinity)になる場合、計算が発散している可能性があります。
  • info < 0: 引数に不正な値が渡されたことを意味します。たとえば、行列の次元が負の値であったり、不正なオプション文字(UPLO, TRANS, DIAGなど)が指定されたりした場合です。
    • トラブルシューティング:
      • trrfs!()に渡す引数の型、次元、値がLAPACKの期待するものと合致しているか、LAPACKの公式ドキュメント(DTRGFS, STRGFS, ZTRGFS, CTRGFSなど、使用しているデータ型に対応するルーチン)を確認してください。
      • 特に、文字列引数(UPLO, TRANS, DIAG)は正しい単一文字である必要があります。
      • 数値引数(行列のサイズNなど)が適切であるか確認してください。

MethodError: no method matching ...

これは、trrfs!()に渡された引数の型や数が、関数が期待するものと一致しない場合に発生します。

  • トラブルシューティング:
    • 引数の型: 例えば、Matrix{Float64}が必要なところにMatrix{Int}を渡したり、Vector{Float64}が必要なところにスカラーを渡したりしていないか確認してください。LAPACK関数は通常、浮動小数点数(Float32, Float64)または複素数(ComplexF32, ComplexF64)の行列とベクトルを期待します。
    • 引数の順序と数: trrfs!()の完全なシグネチャをJuliaのドキュメントまたはソースコードで確認し、全ての必須引数が正しい順序で渡されているかを確認してください。LAPACKの低レベル関数は多くの引数を持ち、一部は作業用配列などになります。
    • 破壊的変更の引数: !が付く関数なので、引数として渡される配列は通常、破壊的に変更されます。配列が変更可能(mutable)な型であるか確認してください(例: Vector, Matrix)。
    • 参照渡し: LAPACKルーチンは、引数として渡された配列のデータを直接操作します。これはJuliaの「参照渡し」に近い動作なので、引数をコピーする必要がある場合は、明示的にコピーを作成してください。

BoundsError: attempt to access ...

これは、配列のインデックスが範囲外になった場合に発生します。LAPACK関数の内部で発生することは稀ですが、作業用配列のサイズが不足している場合や、ユーザーが関数の結果を受け取るために渡した配列のサイズが小さすぎる場合に起こりえます。

  • トラブルシューティング:
    • trrfsに渡す作業用配列(WORK, IWORK)のサイズが、LAPACKドキュメントで推奨されている最小サイズを満たしているか確認してください。これらのサイズは通常、行列の次元に依存します。LAPACKは、最適なパフォーマンスを得るために大きな作業用配列を推奨することがありますが、最小限の作業用配列でも動作は可能です。
    • 出力結果を受け取るための配列(FERR, BERRなど)が、適切な次元とサイズを持っていることを確認してください。

パフォーマンスの問題

直接的なエラーではありませんが、trrfsは反復改良を行うため、非常に悪条件な問題では計算時間が長くなることがあります。

  • トラブルシューティング:
    • 計算時間が異常に長い場合、上記のinfo > 0のエラーと同様に、問題の条件数を確認してください。
    • より安定した、または効率的な前処理(preconditioning)手法を適用できないか検討してください。
    • trrfsは解の精度を「改善」するためのものです。もしも初期の解が非常に悪ければ、いくら反復しても効果が薄い可能性があります。初期解の計算方法を見直すことも考えられます。
  • デバッグ: Juliaのデバッガ(例えば、VS CodeのJulia拡張機能に付属するもの)を使用して、関数呼び出し前後の変数の値を確認し、予期せぬ変更がないか、NaNInfが出現していないかなどをチェックしてください。
  • シンプルな例で試す: まずは、既知の正確な解を持つ非常に小さな行列とベクトルを用いてtrrfs!()を呼び出し、期待通りの動作をするか確認してください。
  • ドキュメントの参照: 最も重要なのは、使用しているJuliaのバージョンに対応するLinearAlgebra.LAPACK.trrfs!()のドキュメント、そして元となるLAPACKのtrrfsルーチン(DTRGFS, STRGFSなど)の公式ドキュメントを詳細に読むことです。引数の意味、期待される型、許容される値、出力の意味などが詳細に記述されています。


一般的なJuliaユーザーが直接この関数を呼び出すことは稀で、通常はより高レベルな関数(例: \演算子やldiv!など)を使用します。しかし、特定の性能要件やデバッグの目的で低レベルなLAPACK関数を直接呼び出す必要がある場合に、trrfs!()は役立ちます。

以下に、LinearAlgebra.LAPACK.trrfs!()の具体的な使用例と、その引数の意味を説明します。

LinearAlgebra.LAPACK.trrfs!()の例

この関数は多くの引数を取るため、非常に複雑に見えます。しかし、それぞれの引数には明確な役割があります。

using LinearAlgebra

# 例として使用する三角形行列 A と右辺ベクトル B
# A * X = B を解くことを想定
N = 3 # 行列のサイズ

# 上三角行列 A を作成
A_upper = [
    3.0  1.0  2.0;
    0.0  4.0  1.0;
    0.0  0.0  5.0
]
A = UpperTriangular(A_upper) # JuliaのUpperTriangular型に変換

# 右辺ベクトル B
B = [10.0; 12.0; 15.0]

# 既に計算されている解 X (ここでは、A \ B で正確な解を求める)
X_exact = A \ B
println("正確な解 X_exact: ", X_exact)

# ここから trrfs! のための準備
# わずかな誤差を含む解 X を作成(trrfs!が改善を試みる)
X_approx = X_exact .+ 0.01 * randn(N) # 少し誤差を加える
println("初期の近似解 X_approx: ", X_approx)

# trrfs! の引数を準備
# これらの引数はLAPACKのDTRRFS/ZTRRFSルーチンの引数に対応します。

# UPLO: 'U' for upper triangular, 'L' for lower triangular
uplo = 'U'

# TRANS: 'N' for A*X=B, 'T' for A'*X=B, 'C' for A^H*X=B
# ここでは A*X=B なので 'N'
trans = 'N'

# DIAG: 'N' for non-unit diagonal, 'U' for unit diagonal
# A の対角成分は1ではないので 'N'
diag = 'N'

# N: size of the matrix A
n = N

# NRHS: number of right-hand sides (number of columns in B and X)
# B がベクトルなので NRHS は 1
nrhs = 1

# A: the matrix A
# LAPACK関数は通常、元の行列の変更された形(factorized form)を受け取るか、
# あるいは元の行列がそのまま使用されます。
# trrfs! の場合、Aは三角行列A自身を保持します。
# JuliaのLAPACKバインディングは、UpperTriangular型から内部のMatrixを取り出して渡す場合があります。
# ここではA_upperをそのまま使用しますが、LAPACKのDTRGFS/ZTRGFSの引数Aは元の三角行列Aそのものです。
# (注意:JuliaのLinearAlgebra.LAPACKモジュールは内部的に適切な変換を行います)

# LD_A: leading dimension of A (stride between columns)
# JuliaのMatrixは列優先なので、LD_Aは行数 N となります
ld_a = N

# B_mat: right-hand side matrix B. Needs to be a matrix for LAPACK.
# ベクトル B を N x 1 の行列に変換
B_mat = reshape(B, N, 1)
ld_b = N

# X_mat: solution matrix X. Needs to be a matrix for LAPACK.
# ベクトル X を N x 1 の行列に変換
X_mat = reshape(X_approx, N, 1)
ld_x = N

# FERR: Output array for forward error bounds. Size (NRHS,).
# NRHS (ここでは1) 個の要素を持つベクトル
ferr = zeros(Float64, nrhs)

# BERR: Output array for backward error bounds. Size (NRHS,).
# NRHS (ここでは1) 個の要素を持つベクトル
berr = zeros(Float64, nrhs)

# WORK: Real workspace array. Size at least 3*N (for DTRGFS).
# DTRGFSの場合、作業領域は3*N以上必要
work = zeros(Float64, 3 * N)

# IWORK: Integer workspace array. Size at least N (for DTRGFS).
iwork = zeros(BlasInt, N) # BlasIntはLAPACKの内部整数型

# trrfs! の呼び出し
# 引数の順番はLAPACKのDTRGFS/ZTRGFSに厳密に従います。
# Juliaのソースコードで LinearAlgebra/src/lapack.jl を確認すると正確なシグネチャが見つかります。
# 例: trrfs!(uplo, trans, diag, A, B, X, ferr, berr, work, iwork)
# 実際には、A, B, X はMatrix{Float64}またはMatrix{ComplexF64}である必要があります。

println("\n--- trrfs!() 呼び出し前 ---")
println("X_mat (変更前): ", X_mat)
println("ferr (変更前): ", ferr)
println("berr (変更前): ", berr)

# LAPACK.trrfs!を呼び出す
# この関数はX_matを破壊的に変更(改良された解で上書き)する可能性があります。
# また、ferrとberrに誤差推定値を書き込みます。
info = LinearAlgebra.LAPACK.trrfs!(uplo, trans, diag, A_upper, B_mat, X_mat, ferr, berr, work, iwork)

println("\n--- trrfs!() 呼び出し後 ---")
println("LAPACK info: ", info) # 0なら成功

if info == 0
    println("改良後の解 X_mat: ", X_mat)
    println("前方誤差 FERR: ", ferr)
    println("後方誤差 BERR: ", berr)

    # 誤差の比較
    println("\n--- 誤差の比較 ---")
    println("初期の解の真の誤差 (||X_approx - X_exact||_inf): ", norm(X_approx - X_exact, Inf))
    println("改良後の解の真の誤差 (||X_mat - X_exact||_inf): ", norm(X_mat - X_exact, Inf))
    println("FERR (推定前方誤差) は真の誤差を上回るべきです。")
else
    println("trrfs! の実行中にエラーが発生しました。info = ", info)
end

コードの説明

  1. using LinearAlgebra: 線形代数モジュールをインポートします。
  2. 行列とベクトルの準備:
    • A_upper: 3x3の上三角行列を通常のMatrix{Float64}として定義します。
    • A = UpperTriangular(A_upper): JuliaのUpperTriangular型に変換することで、線形代数演算が三角行列の構造を考慮して最適化されます。trrfs!には内部の生データであるA_upperを渡します。
    • B: 右辺ベクトル。
    • X_exact: 比較のために、正確な解を\演算子で求めます。
    • X_approx: trrfs!が精度を改善する対象となる、意図的に誤差を含ませた近似解を作成します。
  3. LAPACK引数の設定:
    • uplo, trans, diag: 三角行列の特性(上三角/下三角、転置、対角が単位行列か)を指定する文字コード。これらはLAPACKの標準的な慣習です。
    • n: 行列の次数。
    • nrhs: 右辺の列数。今回はベクトルなので1。
    • A_upper, B_mat, X_mat: LAPACK関数は通常、Fortranベースであり、配列は連続したメモリブロックとして渡されることを期待します。JuliaのMatrixは列優先でメモリに格納されるため、そのまま渡せます。ベクトルはreshapeして1列の行列として渡します。
    • ld_a, ld_b, ld_x: "Leading Dimension"。Juliaの列優先配列では、これは単に行数(N)と同じになります。これは、LAPACKがFortranの配列の次元情報を受け取るために必要です。
    • ferr, berr: 出力用の配列。それぞれ前方誤差と後方残差の推定値が格納されます。nrhs(ここでは1)の要素を持つVector{Float64}として初期化します。
    • work, iwork: LAPACKルーチンが内部計算に使用する一時的な作業用配列。必要な最小サイズはLAPACKのドキュメントで指定されており、通常は行列のサイズに依存します。DTRGFSの場合は3*NNが必要です。
    • BlasInt: LAPACK/BLASライブラリが内部で使用する整数型。JuliaではLinearAlgebra.BlasIntとして提供されます。
  4. LinearAlgebra.LAPACK.trrfs!()の呼び出し:
    • info = LinearAlgebra.LAPACK.trrfs!(...) のように呼び出します。
    • infoは成功(0)またはエラーコードを返します。
  5. 結果の確認:
    • info == 0であれば成功です。
    • X_matは改良された解で上書きされます。
    • ferrberrには誤差推定値が格納されます。ferrは通常、真の誤差の上限を与えるため、norm(X_mat - X_exact, Inf)よりも大きい値になるはずです。

実行結果の解釈

このコードを実行すると、X_approxに加えた小さな誤差がtrrfs!()によってどの程度改善されたか、そしてその誤差がどれくらいだと推定されたかを確認できます。

  • BERR(Backward Error Bound):計算された解 X が、元のシステム AX=B をどの程度正確に満足するかに対応します。つまり、解 X が元のシステム AX′=B′ で真の解 X′ となるような、相対的に小さな摂動 B′ と A′ を持つシステムをどれだけうまく解いたかを示します。厳密には minepsilonmid(A+deltaA)X=(B+deltaB),∣deltaA∣leepsilon∣A∣,∣deltaB∣leepsilon∣B∣ の上限です。
  • FERR(Forward Error Bound):計算された解 X が真の解 X_true からどれくらい離れているかの上限を示します。厳密には frac∣X−X_true∣∗infty∣X∗true∣_infty の上限です。

これらの誤差推定値は、計算された解の信頼性を評価する上で非常に重要です。

  • 高レベル関数との使い分け: ほとんどの場合、JuliaのA \ Bldiv!(Y, A, B)のような高レベルな関数が、内部で最適なLAPACKルーチン(通常はtrtrsのような直接解法)を選択し、必要に応じてリファイメントも行います。trrfs!を直接使うのは、特定の誤差解析が必要な場合や、デバッグ、またはJuliaの標準ライブラリが提供しない特定のLAPACK機能を使いたい場合に限られます。
  • 作業用配列: workiworkは適切なサイズで用意する必要があります。LAPACKのドキュメントを参照してください。
  • 破壊的変更: !が示すように、X_matferrberrは関数内で変更されます。
  • 引数の正確性: LAPACK関数は引数の型、順序、内容に非常に厳格です。少しでも間違えるとLAPACKExceptionMethodErrorが発生します。


\ 演算子(バックスラッシュ演算子)

最も一般的で推奨される方法は、Juliaの組み込みの線形方程式ソルバーである \ 演算子を使用することです。これは、行列の特性(密、疎、対称、正定値など)に基づいて、最適な解法アルゴリズムを自動的に選択します。三角形行列の場合、効率的な前進代入または後退代入が自動的に適用されます。

特徴

  • 要素型に依存: 引数の要素型(Float64, ComplexF64など)に応じて、対応するLAPACKルーチン(DTRTRS, ZTRTRSなど)が呼び出されます。
  • 精度と安定性: 通常、高精度な解を提供し、必要に応じて内部で反復改良やスケーリングを行うことがあります。
  • 自動最適化: 行列の型や特性に応じて、Juliaは最適なLAPACKまたはBLASルーチンを自動的に選択します。
  • 使いやすさ: 最も直感的でシンプルな構文です。


using LinearAlgebra

# 上三角行列 A と右辺ベクトル B
A_upper = [
    3.0  1.0  2.0;
    0.0  4.0  1.0;
    0.0  0.0  5.0
]
A = UpperTriangular(A_upper) # 三角行列型に変換

B = [10.0; 12.0; 15.0]

# 線形方程式 A*X = B を解く
X = A \ B
println("`\` 演算子による解 X: ", X)

# 解の検証(後方残差の計算)
residual = B - A * X
println("残差 (B - A*X): ", residual)
println("残差のノルム: ", norm(residual))

ldiv!() 関数

ldiv!() (left division in-place) 関数は、\ 演算子と同様に線形方程式を解きますが、結果を既存の変数に破壊的に書き込むことで、メモリ割り当てを避けてパフォーマンスを向上させます。特に大規模なシステムやループ内で繰り返し解く場合に有用です。

特徴

  • 因数分解された行列との併用: lu(), cholesky(), qr()などで因数分解された行列オブジェクトと組み合わせて使用できます。
  • 性能: メモリの効率が向上し、多くの場合、\ 演算子よりも高速です。
  • インプレース操作: 解を既存のベクトル/行列に直接書き込むため、新しいメモリ割り当てが不要になります。


using LinearAlgebra

A_upper = [
    3.0  1.0  2.0;
    0.0  4.0  1.0;
    0.0  0.0  5.0
]
A = UpperTriangular(A_upper)

B = [10.0; 12.0; 15.0]

# 解を格納するベクトルを初期化 (Bと同じ型とサイズ)
X_ldiv = similar(B)

# ldiv! で解を計算
ldiv!(X_ldiv, A, B)
println("`ldiv!()` による解 X_ldiv: ", X_ldiv)

因数分解を用いた解法 (lu(), cholesky(), qr() など)

行列を事前に因数分解( factorization )しておくと、同じ行列 A で異なる右辺ベクトル B の線形方程式を繰り返し解く場合に非常に効率的です。因数分解オブジェクトは、\ 演算子や ldiv!() と組み合わせて使用できます。

特徴

  • タイプに応じた選択: 三角行列に対しては、その構造を活かした因数分解(例: lu()は通常、一般の密行列に適用されますが、三角形行列の場合も効率的に機能します)が適用されます。
  • 安定性: 因数分解の過程で、行列の特異性や条件数に関する情報が得られる場合があります。
  • 効率性: 因数分解は行列の構造に特化した最適化されたアルゴリズムを使用します。一度因数分解すれば、複数の右辺ベクトルに対して高速に解を求められます。


using LinearAlgebra

A_upper = [
    3.0  1.0  2.0;
    0.0  4.0  1.0;
    0.0  0.0  5.0
]
A = UpperTriangular(A_upper)

B1 = [10.0; 12.0; 15.0]
B2 = [1.0; 2.0; 3.0]

# LU分解(三角形行列なので、UはAとほとんど同じになります)
# 一般的なLU分解は、UpperTriangular型をLUFactorization型に変換します。
# 三角行列なので、内部ではtrtrsのような解法が使われます。
F = lu(A)

# 因数分解オブジェクト F を使って解く
X1 = F \ B1
X2 = F \ B2

println("LU分解と`\` 演算子による解 X1: ", X1)
println("LU分解と`\` 演算子による解 X2: ", X2)

# ldiv! と組み合わせる
X_ldiv_B1 = similar(B1)
ldiv!(X_ldiv_B1, F, B1)
println("LU分解と`ldiv!()` による解 X_ldiv_B1: ", X_ldiv_B1)

trrfs!() が提供する前方誤差 (FERR) や後方残差 (BERR) の推定値は、高レベルな関数では直接提供されません。しかし、解の品質を評価するための代替手段は存在します。

  • 多倍長精度計算: 非常に悪条件な問題で、通常のFloat64の精度では不十分な場合、BigFloatなどの多倍長精度を使用することで、より正確な解を得ることができます。これは計算コストが高くなりますが、高い精度が絶対的に必要な場合に有効です。

    A_big = BigFloat[
        3.0  1.0  2.0;
        0.0  4.0  1.0;
        0.0  0.0  5.0
    ]
    B_big = BigFloat[10.0; 12.0; 15.0]
    
    X_big = UpperTriangular(A_big) \ B_big
    println("BigFloatによる解 X_big: ", X_big)
    
  • 条件数 (cond() ): 行列の条件数 $ \kappa(A) = |A| |A^{-1}| $ は、線形方程式 $ AX=B $ の問題がどの程度「悪条件」であるかを示す指標です。条件数が大きいほど、入力の小さな摂動(浮動小数点誤差も含む)が解に大きな誤差を引き起こす可能性があります。条件数が大きい場合、trrfs!()のような反復改良が役立つことがありますが、そうでなければ必要ないかもしれません。

    println("行列 A の無限大ノルム条件数: ", cond(A, Inf))
    println("行列 A の2ノルム条件数: ", cond(A, 2))
    
  • 残差ノルム: 最も直接的な方法は、計算された解 X を使って残差 R=B−AX を計算し、そのノルム ∣R∣ を調べることです。このノルムが小さければ小さいほど、解は正確です(後方誤差が小さい)。

    X = A \ B
    residual = B - A * X
    println("残差のノルム: ", norm(residual))
    

ほとんどの場合、LinearAlgebra.LAPACK.trrfs!()を直接使用する代わりに、Juliaの標準ライブラリが提供する高レベルな機能(\ 演算子、ldiv!()、因数分解)を使用することを強くお勧めします。これらはより安全で、使いやすく、多くの場合、内部でLAPACKの最適化されたルーチン(必要に応じて誤差改良も含む)を活用して、十分な性能と精度を提供します。