Please note that the recommended version of Scilab is 2025.0.0. This page might be outdated.
However, this page did not exist in the previous stable version.
pcg
プリコンディショナ付き共役勾配法
呼び出し手順
[x, flag, err, iter, res] = pcg(A, b [, tol [, maxIter [, M [, M2 [, x0 [, verbose]]]]]]) [x, flag, err, iter, res] = pcg(A, b [key=value,...])
パラメータ
- A
行列, または指定した
x
についてA*x
を計算する関数またはリスト. 以下にAの型に応じたA*xの計算に関する説明を示します.行列.
Aが行列の場合, 通常の行列 または疎行列とすることができます関数.
Aが関数の場合, 以下のヘッダーを有する必要があります :function y=A(x)
list.
Aがリストの場合, リストの最初の要素は関数で,リストのその他の 要素(インデックス2から末尾まで)は関数の引数となります. 関数がコールされる時, xのカレントの値は関数に最初の引数として渡されます. 他の引数はリストに指定された値です.
- b
右辺ベクトルr (大きさ: nx1)
- tol
相対許容誤差 (デフォルト: 1e-8). 終了条件は残差r=b-Axの2次ノルムを 右辺bの2次ノルムで割った値に基づきます.
- maxIter
反復最大回数 (デフォルト: n)
- M
プリコンディショナ: 通常の行列または疎行列または
M\x
を返す関数 (デフォルト: none)- M2
プリコンディショナ: 通常の行列または疎行列または
x
についてM2\x
を返す関数 (デフォルト: none)- x0
初期推定ベクトル (デフォルト: zeros(n,1))
- verbose
冗長なログを有効にする場合は1に設定(デフォルト 0)
- x
解ベクトル
- flag
maxi
回の反復回数内にpcg
が指定した 許容誤差に収束する場合に 0, そうでない場合に 1- err
最終残差相対ノルム (右辺bの2次ノルムが使用されます)
- iter
実行された反復回数
- res
残差相対ノルムのベクトル
説明
プリコンディショナ有りまたは無しの
共役勾配法により線形システム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)
Aを疎行列または関数またはリストとして指定する例
以下の例では,疎行列を同時に処理する方法を示します. 右辺を計算する関数が "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構文の例
以下の例ではkey=value"構文により引数を指定する方法を示します. これにより,位置を固定せずに引数を設定できるようになり, 引数のリストにおいてその順番に依存せずに引数を設定できます. 利用可能なキーはオプションの引数の名前で,以下に示します: tol, maxIter, %M, %M2, x0, verbose. 以下の例では, maxIterオプションの前にverboseオプションを指定します. "key=value" 構文ではなく, 位置を指定する引数の場合は, maxIterを最初にverboseを後に指定する必要があります.
// 引数をkey=value構文で指定する例 A=[100,1;1,10]; b=[101;11]; [xcomputed, flag, err, iter, res]=pcg(A,b,verbose=1); // key=value構文では, 順番は関係ありません [xcomputed, flag, err, iter, res]=pcg(A,b,verbose=1,maxIter=0);
参照
作者
Sage Group, IRISA, 2004
Serge Steer, INRIA, 2006
Michael Baudin, INRIA, 2008-2009
参考文献
"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
<< nnz | Sparses Matrix | qmr >> |