プリコンディショナ付き共役勾配法
[x, flag, err, iter, res] = pcg(A, b [, tol [, maxIter [, M [, M2 [, x0 [, verbose]]]]]]) [x, flag, err, iter, res] = pcg(A, b [key=value,...])
行列, または指定したx
について
A*x
を計算する関数またはリスト.
以下にAの型に応じたA*xの計算に関する説明を示します.
行列.
Aが行列の場合, 通常の行列
または疎行列とすることができます
関数.
Aが関数の場合,
以下のヘッダーを有する必要があります :
list.
Aがリストの場合,
リストの最初の要素は関数で,リストのその他の
要素(インデックス2から末尾まで)は関数の引数となります.
関数がコールされる時,
xのカレントの値は関数に最初の引数として渡されます.
他の引数はリストに指定された値です.
右辺ベクトルr (大きさ: nx1)
相対許容誤差 (デフォルト: 1e-8). 終了条件は残差r=b-Axの2次ノルムを 右辺bの2次ノルムで割った値に基づきます.
反復最大回数 (デフォルト: n)
プリコンディショナ: 通常の行列または疎行列または
M\x
を返す関数 (デフォルト: none)
プリコンディショナ: 通常の行列または疎行列または
x
についてM2\x
を返す関数 (デフォルト:
none)
初期推定ベクトル (デフォルト: zeros(n,1))
冗長なログを有効にする場合は1に設定(デフォルト 0)
解ベクトル
maxi
回の反復回数内に
pcg
が指定した
許容誤差に収束する場合に 0, そうでない場合に 1
最終残差相対ノルム (右辺bの2次ノルムが使用されます)
実行された反復回数
残差相対ノルムのベクトル
プリコンディショナ有りまたは無しの
共役勾配法により線形システムAx=b
を解きます.
プリコンディショナは,
対称正定行列M
またはM=M1*M2
となるような
2つの行列M1
およびM2
により定義されます.
この場合,関数はinv(M)*A*x = inv(M)*b
を
x
について解きます.
M
, M1
および
M2
は,対応する左除算
y=Mi\x
を計算する
呼び出し手順y=Milx(x)
を有する
Scilab関数とすることができます.
A
行列は対称正定行列
(通常の行列または疎行列)または
y=A*x
を計算する
呼び出し手順y=Ax(x)
を有する関数と
する必要があります.
以下の例では,2つの線形システムを解きます. 最初の行列の条件数は ~0.02 に等しくなり, アルゴリズムはちょうど10回の反復で収束します. これが行列の大きさの場合,共役勾配法で指定される動作となります. 後者は,条件数1.d-6に等しくなり,22回の反復で収束します. これは,パラメータ maxIterを30に設定する理由です. "key=value" 構文のその他の例については以下を参照ください.
//良い条件の例 A=[ 94 0 0 0 0 28 0 0 32 0 0 59 13 5 0 0 0 10 0 0 0 13 72 34 2 0 0 0 0 65 0 5 34 114 0 0 0 0 0 55 0 0 2 0 70 0 28 32 12 0 28 0 0 0 0 87 20 0 33 0 0 0 0 0 28 20 71 39 0 0 0 10 0 0 32 0 39 46 8 0 32 0 0 0 12 33 0 8 82 11 0 0 65 55 0 0 0 0 11 100]; b=ones(10,1); [x, fail, err, iter, res]=pcg(A,b,1d-12,15); mprintf(" fail=%d, iter=%d, errrel=%e\n",fail,iter,err) //悪い条件の例 A=[ 894 0 0 0 0 28 0 0 1000 70000 0 5 13 5 0 0 0 0 0 0 0 13 72 34 0 0 0 0 0 6500 0 5 34 1 0 0 0 0 0 55 0 0 0 0 70 0 28 32 12 0 28 0 0 0 0 87 20 0 33 0 0 0 0 0 28 20 71 39 0 0 0 0 0 0 32 0 39 46 8 0 1000 0 0 0 12 33 0 8 82 11 70000 0 6500 55 0 0 0 0 11 100]; [x, fail, err, iter, res]=pcg(A,b,maxIter=30,tol=1d-12); mprintf(" fail=%d, iter=%d, errrel=%e\n",fail,iter,err) | ![]() | ![]() |
以下の例では,疎行列を同時に処理する方法を示します. 右辺を計算する関数が "pcg" プリミティブに指定される 場合も示します. この例における最後のケースでは, リストがプリミティブに指定される場合です.
//良い条件の例 A=[ 94 0 0 0 0 28 0 0 32 0 0 59 13 5 0 0 0 10 0 0 0 13 72 34 2 0 0 0 0 65 0 5 34 114 0 0 0 0 0 55 0 0 2 0 70 0 28 32 12 0 28 0 0 0 0 87 20 0 33 0 0 0 0 0 28 20 71 39 0 0 0 10 0 0 32 0 39 46 8 0 32 0 0 0 12 33 0 8 82 11 0 0 65 55 0 0 0 0 11 100]; b=ones(10,1); // Aを疎行列に変換 Asparse=sparse(A); [x, fail, err, iter, res]=pcg(Asparse,b,maxIter=30,tol=1d-12); mprintf(" fail=%d, iter=%d, errrel=%e\n",fail,iter,err) // 右辺を計算する関数を定義. function y=Atimesx(x) A=[ 94 0 0 0 0 28 0 0 32 0 0 59 13 5 0 0 0 10 0 0 0 13 72 34 2 0 0 0 0 65 0 5 34 114 0 0 0 0 0 55 0 0 2 0 70 0 28 32 12 0 28 0 0 0 0 87 20 0 33 0 0 0 0 0 28 20 71 39 0 0 0 10 0 0 32 0 39 46 8 0 32 0 0 0 12 33 0 8 82 11 0 0 65 55 0 0 0 0 11 100]; y=A*x endfunction // スクリプトAtimesx をプリミティブに指定 [x, fail, err, iter, res]=pcg(Atimesx,b,maxIter=30,tol=1d-12); mprintf(" fail=%d, iter=%d, errrel=%e\n",fail,iter,err) // 右辺を計算する関数を定義. function y=Atimesxbis(x, A) y=A*x endfunction // リストをプリミティブに指定 Alist = list(Atimesxbis,Asparse); [x, fail, err, iter, res]=pcg(Alist,b,maxIter=30,tol=1d-12); mprintf(" fail=%d, iter=%d, errrel=%e\n",fail,iter,err) | ![]() | ![]() |
以下の例ではkey=value"構文により引数を指定する方法を示します. これにより,位置を固定せずに引数を設定できるようになり, 引数のリストにおいてその順番に依存せずに引数を設定できます. 利用可能なキーはオプションの引数の名前で,以下に示します: tol, maxIter, %M, %M2, x0, verbose. 以下の例では, maxIterオプションの前にverboseオプションを指定します. "key=value" 構文ではなく, 位置を指定する引数の場合は, maxIterを最初にverboseを後に指定する必要があります.
"Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra, Eijkhout, Pozo, Romine, and Van der Vorst, SIAM Publications, 1993, ftp netlib2.cs.utk.edu/linalg/templates.ps
"Iterative Methods for Sparse Linear Systems, Second Edition", Saad, SIAM Publications, 2003, ftp ftp.cs.umn.edu/dept/users/saad/PS/all_ps.zip
Version | Description |
5.5.1 | pcg was removed after Scilab 5.5.1.
conjgrad replaces it. |