Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Please login or create an account
Change language to: English - Français - Português - Русский

Please note that the recommended version of Scilab is 6.0.2. This page might be outdated.
See the recommended documentation of this function

Scilabヘルプ >> Sparses Matrix > Linear Equations (Iterative Solvers) > conjgrad

conjgrad

共役勾配ソルバ

呼び出し手順

[x, flag, err, iter, res] = conjgrad(A, b [, method [, tol [, maxIter [, M [, M2 [, x0 [, verbose]]]]]]])
[x, flag, err, iter, res] = conjgrad(A, b [, method [, key=value,...]])

引数

A

指令した各xについてA*xを計算する 行列または関数またはリスト. Aのそれぞれの型に関する A*x の計算に関して以下に示します.

  • 行列.Aが行列の場合, 通常または疎行列が使用可能

  • 関数.Aが関数の場合, 以下のヘッダを有する 必要があります :

    function y=A(x)
  • リスト.Aがリストの場合, リストの最初の要素には関数を指定し, 添字2から最後までのリストの他の要素には関数の引数を指定します. この関数がコールされた際, xのカレントの値が関数の最初の引数に指定され, その他の引数にはリストで指定されたものが指定されます.

b

右辺ベクトル (大きさ: nx1)

mehtod

スカラー文字列, "pcg", "cgs", "bicg" または "bicgstab" (デフォルト "bicgstab")

tol

相対許容誤差 (デフォルト: 1e-8). 終了基準は残差 r=b-Ax の二乗ノルムを右辺 b の二乗ノルムで割ったものに 基づきます.

maxIter

最大反復回数 (デフォルト: n)

M

プリコンディショナ: 通常または疎行列または M\x を返す関数 (デフォルト: none)

M2

プリコンディショナ: 通常または疎行列または各xM2\x を返す関数 (デフォルト: none)

x0

初期推定ベクトル (デフォルト: zeros(n,1))

verbose

冗長なログを有効にする場合は1を指定 (デフォルト 0)

x

解ベクトル

flag

conjgrad が反復回数maxi以内に 許容誤差以内に収束した場合は 0, それ以外は 1

err

残差ノルムの最終値 (右辺 b の二乗ノルムを使用)

iter

実行した反復回数

res

相対ノルム残差のベクトル

説明

プリコンディショニング有りまたは無しの 共役勾配法により線形システムAx=b を解きます. プリコンディショナは対称正定行列M,または M=M1*M2となる2つの行列 M1M2 により定義されます. この場合,この関数はinv(M)*A*x = inv(M)*bxについて解きます. M, M1 および M2 は, 呼び出し手順y=Milx(x) で対応する左除算y=Mi\xを計算するScilab関数と することができます.

入力引数 method は,使用するソルバを指定します:

  • "pcg" プリコンディショナ付き共役勾配法
  • "cgs" プリコンディショナ付き二乗共役勾配法
  • "bicg" プリコンディショナ付き双共役勾配法
  • "bicgstab" プリコンディショナ付き安定化双共役勾配法 (デフォルト)

method="pcg"の場合, A 行列は 対称正定行列(通常または疎行列)または呼び出し手順y=Ax(x)y=A*xを計算する関数とする必要があります. その他の場合 (method="cgs", "bicg" or "bicgstab"), A は正方行列であることのみ必要です.

良条件および悪条件の問題の例

以下の例では, 2つの線形システムを解きます. 最初の行列は条件数が ~0.02で,アルゴリズムは10回で収束します. これは行列の大きさと同じで,共役勾配法で期待される動作です. 2番目の行列は条件数が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] = conjgrad(A, b, "bicg", 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] = conjgrad(A, b, method="pcg", maxIter=30, tol=1d-12);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)

Aに疎行列または関数またはリストを指定する例

以下の例は,この手法が疎行列を同様に処理できることを示すものです. また,この例は右辺を計算する関数を"conjgrad"のプリミティブに指定する ケースも示します. この例の最後のケースは, リストをプリミティブに指定すた場合です.

// 良条件問題
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] = conjgrad(Asparse, b, "cgs", 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] = conjgrad(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] = conjgrad(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 です. 以下の例では, verbose オプションが maxIterオプションの 前に指定されていることに注意してください. "key=value"構文ではない場合は引数の位置が固定されるため, maxIter は先, verbose は後で使用する必要があります.

// key=value 構文で指定される引数の例
A = [100 1; 1 10];
b = [101; 11];
[xcomputed, flag, err, iter, res] = conjgrad(A, b, verbose=1);
// key=value 構文の場合, 順序は関係ありません
[xcomputed, flag, err, iter, res] = conjgrad(A, b, verbose=1, maxIter=0);

参照

  • backslash — (\) 左行列除算.
  • qmr — プリコンディショナ付きのQuasi Minimal Residual法
  • gmres — Generalized Minimum RESidual 法

参考文献

PCG

"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

CGS

"CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems" by Peter Sonneveld.

Original article

Article on ACM

Some theory around CGS

BICG

"Numerical Recipes: The Art of Scientific Computing." (third ed.) by William Press, Saul Teukolsky, William Vetterling, Brian Flannery.

http://apps.nrbook.com/empanel/index.html?pg=87

Article on ACM

Some theory around BICG

BICGSTAB

"Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems" by Henk van der Vorst. 339

Original article

Article on ACM

Some theory around BICG

履歴

VersionDescription
5.5.0 導入
Scilab Enterprises
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors
Last updated:
Mon Feb 12 23:12:40 CET 2018