◆ Dpotrf()

Sub Dpotrf	(	Uplo As	String,
		N As	Long,
		A() As	Double,
		Info As	Long
	)

係数行列のコレスキー分解 (正定値対称行列)

目的: 本ルーチンは正定値対称行列Aのコレスキー分解を計算する. 分解は次の形式である.
A = U^T * U (Uplo = "U"の場合)

A = L * L^T (Uplo = "L"の場合)

ここで, Uは上三角行列, Lは下三角行列である.
本ルーチンはLevel 3 BLASを呼び出すブロック版のアルゴリズムを使用する.

引数

[in]	Uplo	= "U": Aの上三角部分を格納. = "L": Aの下三角部分を格納.
[in]	N	行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in,out]	A()	配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] N×N正定値対称行列 A. Uploに従い上三角部分あるいは下三角部分が参照される. [out] Info = 0の場合, コレスキー分解 A = U^TU または A = LL^T のUまたはL.
[out]	Info	= 0: 正常終了. = -1: パラメータ Uplo の誤り. (Uplo != "U"および"L") = -2: パラメータ N の誤り. (N < 0) = -3: パラメータ A() の誤り. = i > 0: i×i首座小行列が正定値でないため分解を完了できなかった.

出典: LAPACK

使用例: 連立一次方程式 Ax = B を解き, 同時にAの条件数の逆数の推定値(RCond)を求める.
ただし, Aは正定値対称行列で
( 2.2 -0.11 -0.32 ) ( -1.566 )

A = ( -0.11 2.93 0.81 ), B = ( -2.8425 )

( -0.32 0.81 2.37 ) ( -1.1765 )

とする.
Sub Ex_Dpotrf()

Const N = 3

Dim A(N - 1, N - 1) As Double, B(N - 1) As Double

Dim ANorm As Double, RCond As Double, Info As Long

A(0, 0) = 2.2

A(1, 0) = -0.11: A(1, 1) = 2.93

A(2, 0) = -0.32: A(2, 1) = 0.81: A(2, 2) = 2.37:

B(0) = -1.566: B(1) = -2.8425: B(2) = -1.1765

ANorm = Dlansy("1", "L", N, A())

Call Dpotrf("L", N, A(), Info)

If Info = 0 Then Call Dpotrs("L", N, A(), B(), Info)

If Info = 0 Then Call Dpocon("L", N, A(), ANorm, RCond, Info)

Debug.Print "X =", B(0), B(1), B(2)

Debug.Print "RCond =", RCond

Debug.Print "Info =", Info

End Sub

Dlansy
Function Dlansy(Norm As String, Byval Uplo As String, N As Long, A() As Double, Optional Info As Long) As Double
行列の1ノルム, フロベニウスノルム, 無限ノルム, または, 要素の最大絶対値 (対称行列)

Dpocon
Sub Dpocon(Uplo As String, N As Long, A() As Double, ANorm As Double, RCond As Double, Info As Long)
行列の条件数 (正定値対称行列)

Dpotrs
Sub Dpotrs(Uplo As String, N As Long, A() As Double, B() As Double, Info As Long, Optional Nrhs As Long=1)
分解済の連立一次方程式 AX = B の解 (正定値対称行列)

Dpotrf
Sub Dpotrf(Uplo As String, N As Long, A() As Double, Info As Long)
係数行列のコレスキー分解 (正定値対称行列)

実行結果: X = -0.8 -0.92 -0.29

RCond = 0.446791078068956

Info = 0