XLPack 7.0
XLPack 数値計算ライブラリ (Excel VBA) リファレンスマニュアル
読み取り中…
検索中…
一致する文字列を見つけられません

◆ Dgbsv()

Sub Dgbsv ( N As  Long,
Kl As  Long,
Ku As  Long,
Ab() As  Double,
IPiv() As  Long,
B() As  Double,
Info As  Long,
Optional Nrhs As  Long = 1 
)

(シンプルドライバ) 連立一次方程式 AX = B の解 (一般帯行列)

目的
本ルーチンは次の連立一次方程式を解く.
A * X = B
ここで, Aは下帯幅kl, 上帯幅kuのn×n帯行列, また, XおよびBはn×nrhs行列である.

まず, 行交換によるピボットの部分選択を行うLU分解を用いて, 次のようにAを分解する.
A = L * U
ここで, Lは置換行列と下帯幅klで対角要素が1の下三角行列の積, そして, Uは上帯幅kl+kuの上三角行列である. 次に, 分解されたAを用いて連立方程式 A * X = B の解を求める.
引数
[in]N連立方程式の数, すなわち, 行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]KlAの下帯幅. (Kl >= 0)
[in]KuAの上帯幅. (Ku >= 0)
[in,out]Ab()配列 Ab(LAb1 - 1, LAb2 - 1) (LAb1 >= 2Kl + Ku + 1, LAb2 >= N)
[in] 帯行列形式の行列Aを第Kl+1〜2Kl+Ku+1行に与える. 第1〜Kl行は設定する必要がない.
[out] 分解結果: Uは上帯幅Kl+Kuの上三角行列として第1〜Kl+Ku+1行に格納される. また, 分解中に使われた乗数が第Kl+Ku+2〜2Kl+Ku+1行に格納される. 詳細は下記を参照のこと.
[out]IPiv()配列 IPiv(LIPiv - 1) (LIPiv >= N)
置換行列Pを定義するピボットインデックス. 第i行が第IPiv(i-1)行と交換されたことを表す.
[in,out]B()配列 B(LB1 - 1, LB2 - 1) (LB1 >= max(1, N), LB2 >= Nrhs) (2次元配列) または B(LB - 1) (LB >= max(1, N), Nrhs = 1) (1次元配列)
[in] N×Nrhs右辺行列 B.
[out] Info = 0 の場合, N×Nrhs解行列 X.
[out]Info= 0: 正常終了.
= -1: パラメータ N の誤り. (N < 0)
= -2: パラメータ Kl の誤り. (Kl < 0)
= -3: パラメータ Ku の誤り. (Ku < 0)
= -4: パラメータ Ab() の誤り.
= -5: パラメータ IPiv() の誤り.
= -6: パラメータ B() の誤り.
= -8: パラメータ Nrhs の誤り. (nrhs < 0)
= i > 0: Uのi番目の対角要素が0である. 分解を完了したが, Uが特異であり解は計算されなかった.
[in]Nrhs(省略可)
右辺の数, すなわち, 行列Bの列数. (Nrhs >= 0) (Nrhs = 0 の場合, 処理を行わずに戻る) (省略時 = 1)
詳細
次の例は, N = 6, Kl = 2, Ku = 1 の場合の帯行列形式を表す.
入力:
   *    *    *    +    +    +
   *    *    +    +    +    +
   *   a12  a23  a34  a45  a56
  a11  a22  a33  a44  a55  a66
  a21  a32  a43  a54  a65   *
  a31  a42  a53  a64   *    *

出力:
   *    *    *   u14  u25  u36
   *    *   u13  u24  u35  u46
   *   u12  u23  u34  u45  u56
  u11  u22  u33  u44  u55  u66
  m21  m32  m43  m54  m65   *
  m31  m42  m53  m64   *    *
*で示された配列要素は使用されない. +で示された要素は入力不要であるが, 行交換の結果フィルインが生じるUの要素を格納するために必要である.
出典
LAPACK
使用例
連立一次方程式 Ax = B を解き, 同時にAの条件数の逆数の推定値(RCond)を求める. ただし,
( 2.34 0.57 0 ) ( 0.7416 )
A = ( 0.65 1.98 -1.39 ), B = ( 0.7885 )
( 0 1.50 1.73 ) ( 1.0833 )
とする.
Sub Ex_Dgbsv()
Const N = 3, Kl = 1, Ku = 1
Dim Ab(2 * Kl + Ku, N - 1) As Double, B(N - 1) As Double, IPiv(N - 1) As Long
Dim ANorm As Double, RCond As Double, Info As Long
Ab(1, 1) = 0.57: Ab(1, 2) = -1.39
Ab(2, 0) = 2.34: Ab(2, 1) = 1.98: Ab(2, 2) = 1.73
Ab(3, 0) = 0.65: Ab(3, 1) = 1.5
B(0) = 0.7416: B(1) = 0.7885: B(2) = 1.0833
ANorm = Dlangb("1", N, Kl, Ku, Ab(), , Kl)
Call Dgbsv(N, Kl, Ku, Ab(), IPiv(), B(), Info)
If Info = 0 Then Call Dgbcon("1", N, Kl, Ku, Ab(), IPiv(), ANorm, RCond, Info)
Debug.Print "X =", B(0), B(1), B(2)
Debug.Print "RCond =", RCond
Debug.Print "Info =", Info
End Sub
Function Dlangb(Norm As String, N As Long, Kl As Long, Ku As Long, Ab() As Double, Optional Info As Long, Optional Offset As Long=0) As Double
行列の1ノルム, フロベニウスノルム, 無限ノルム, または, 要素の最大絶対値 (帯行列)
Sub Dgbcon(Norm As String, N As Long, Kl As Long, Ku As Long, Ab() As Double, IPiv() As Long, ANorm As Double, RCond As Double, Info As Long)
行列の条件数 (一般帯行列)
Sub Dgbsv(N As Long, Kl As Long, Ku As Long, Ab() As Double, IPiv() As Long, B() As Double, Info As Long, Optional Nrhs As Long=1)
(シンプルドライバ) 連立一次方程式 AX = B の解 (一般帯行列)
実行結果
X = 0.2 0.48 0.21
RCond = 0.364187455306431
Info = 0