18.5 特殊なソルバー
: x = bicg (A, b)
: x = bicg (A, b, tol)
: x = bicg (A, b, tol, maxit)
: x = bicg (A, b, tol, maxit, M)
: x = bicg (A, b, tol, maxit, M1, M2)
: x = bicg (A, b, tol, maxit, M, [], x0)
: x = bicg (A, b, tol, maxit, M1, M2, x0)
: x = bicg (A, b, tol, maxit, M, [], x0, …)
: x = bicg (A, b, tol, maxit, M1, M2, x0, …)
: [x, flag, relres, iter, resvec] = bicg (A, b, …)
双共役勾配反復法を使用して 線形方程式系を解きます 。A * x = b
入力引数は次のとおりです。
Aは線形システムの行列であり、正方行列でなければなりません。 A は、 および と Afcnなるような行列、関数ハンドル、またはインライン関数として渡すことができます 。 x0の後にへの追加パラメータ を渡すこともできます。 Afcn (x, "notransp") = A * xAfcn (x, "transp") = A' * xAfcn
bは右側のベクトルです。Aと同じ行数の列ベクトルでなければなりません。
tol は残差誤差に必要な相対許容値です 。 ≤ の場合、反復は停止します 。tolが省略されるか空の場合、許容値 1e-6 が使用されます。 b - A * xnorm (b - A * x)tol * norm (b)
maxit は許容される最大反復回数です。maxit が省略されるか空の場合、 値 20 が使用されます。
M1、M2は前処理子です。前処理子M はとして与えられます。M1とM2 は両方とも、 または および または となる ような 行列、または関数ハンドルまたはインライン関数として渡すことができます 。M1が省略されるか空の場合、前処理は適用されません。前処理されたシステムは、理論的には、メソッドを線形システムに 適用し、 を設定すること と同等です 。 M = M1 * M2gg (x, "notransp") = M1 \ xg (x, "notransp") = M2 \ xg (x, "transp") = M1' \ xg (x, "transp") = M2' \ xbicginv (M1) * A * inv (M2) * y = inv (M1) * binv (M2') * A' * inv (M1') * z = inv (M2') * bx = inv (M2) * y
x0は初期推定値です。x0 が省略または空の場合、関数はデフォルトでx0 をゼロ ベクトルに設定します。
x0に続く引数はすべてパラメータとして扱われ、に与えられた関数 ( AfcnまたはMfcnbicg ) または のいずれかに適切な方法で渡されます。
出力パラメータは次のとおりです。
x はの解の計算された近似値です 。アルゴリズムが収束しなかった場合、xは残差が最小となる反復です。 A * x = b
フラグは終了ステータスを示します:
0: アルゴリズムは規定の許容範囲内に収束しました。
1: アルゴリズムが収束せず、反復の最大回数に達しました。
2: 前処理行列は特異である。
3: アルゴリズムが停滞しました。つまり、現在の反復xと前の反復の差の絶対値が 未満です。 eps * norm (x,2)
4: 中間値が小さすぎたり大きすぎたりして信頼性の高い計算が行えなくなったため、アルゴリズムを続行できませんでした。
relres は、ユークリッドノルムで測定された最終残差と初期値の比です。
iter はxが計算される反復です。
resvec は各反復における残差を含むベクトルです。実行される反復の合計数は で与えられます 。 length (resvec) - 1
三角行列に関する簡単な問題を考えてみましょう
n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afcn = @(x, string) strcmp (string, "notransp") * (A * x) + ...
strcmp (string, "transp") * (A' * x);
Mfcn = @(x, string) strcmp (string, "notransp") * (M \ x) + ...
strcmp (string, "transp") * (M' \ x);
M1fcn = @(x, string) strcmp (string, "notransp") * (M1 \ x) + ...
strcmp (string, "transp") * (M1' \ x);
M2fcn = @(x, string) strcmp (string, "notransp") * (M2 \ x) + ...
strcmp (string, "transp") * (M2' \ x);
例1:最も簡単な使用法bicg x = bicg (A, b)
例2: 計算bicgする関数を使用して A*xA'*x x = bicg (Afcn, b, [], n)
例3: bicg前処理行列M x = bicg (A, b, 1e-6, n, M)
例4: bicg関数を前処理として用いる x = bicg (Afcn, b, 1e-6, n, Mfcn)
例5: bicg前処理行列M1 とM2 x = bicg (A, b, 1e-6, n, M1, M2)
例6: bicg関数を前処理として用いる x = bicg (Afcn, b, 1e-6, n, M1fcn, M2fcn)
例7: bicg引数を必要とする関数を入力としてargument
function y = Ap (A, x, string, z)
## compute A^z * x or (A^z)' * x
y = x;
if (strcmp (string, "notransp"))
for i = 1:z
y = A * y;
endfor
elseif (strcmp (string, "transp"))
for i = 1:z
y = A' * y;
endfor
endif
endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p); x = bicg (Apfcn, b, [], [], [], [], [], 2);
参照:
Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM.
See also: bicgstab, cgs, gmres, pcg, qmr, tfqmr.
: x = bicgstab (A, b, tol, maxit, M1, M2, x0, …)
: x = bicgstab (A, b, tol, maxit, M, [], x0, …)
: [x, flag, relres, iter, resvec] = bicgstab (A, b, …)
A x = b安定化された双共役勾配反復法を使用して 解きます。
入力パラメータは次のとおりです。
Aは線形システムの行列であり、正方行列でなければなりません。 A はAfcn、のような行列、関数ハンドル、またはインライン関数として渡すことができますAfcn(x) = A * x。 への追加パラメータは、x0Afcnの後に渡されます。
bは右側のベクトルです。Aと同じ行数の列ベクトルでなければなりません。
tol は残差誤差に必要な相対許容値です 。 ≤ の場合、反復は停止します 。tolが省略されるか空の場合、許容値 1e-6 が使用されます。 b - A * xnorm (b - A * x)tol * norm (b)
maxit は外部反復の最大回数です。指定されないか [] に設定されている場合、デフォルト値min (20, numel (b))が使用されます。
M1、M2は前処理子です。前処理子 Mは として与えられます。M1とM2 は両方とも、または となる ような 行列、関数ハンドル、またはインライン関数として渡すことができます。使用される手法は適切な前処理です。つまり、 を解いてから と なります。 M = M1 * M2gg(x) = M1 \ xg(x) = M2 \ xA * inv (M) * y = bx = inv (M) * y
x0初期推定値。指定されないか [] に設定されている場合、デフォルト値が使用されます。 zeros (size (b))
x0に続く引数はパラメータとして扱われ、に渡される関数 ( AまたはMbicstab ) のいずれかに適切な方法で渡されます。
出力パラメータは次のとおりです。
xは計算された近似値です。この方法が収束しない場合は、最小残差で反復されます。
フラグは終了ステータスを示します:
0: 反復は選択された許容範囲内に収束しました
1: 収束前に最大反復回数に達した
2: 前処理行列は特異である
3: アルゴリズムが停滞した
4: ゼロ除算のためアルゴリズムを続行できません
relres は、 として得られた相対残差です 。 (A*x-b) / norm(b)
iterはxが計算される(おそらく半分の)反復です。半分の反復の場合はiter + 0.5
resvec は、各半分の反復と全体の反復の残差を含むベクトルです (各反復でx が2 つのステップで計算されるため、半分の反復もあります)。これにより、実行された (合計) 反復の合計数を確認できます。 (length(resvec) - 1) / 2
三角行列に関する簡単な問題を考えてみましょう
n = 20; A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... sparse (1, 2, 1, 1, n) * n / 2); b = A * ones (n, 1); restart = 5; [M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A) M = M1 * M2; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(x) M2 \ x;
例1:最も簡単な使用法bicgstab x = bicgstab (A, b, [], n)
例2: bicgstab計算する関数 A * x x = bicgstab (Afcn, b, [], n)
例3: bicgstab前処理行列M x = bicgstab (A, b, [], 1e-06, n, M)
例4: bicgstab関数を前処理として用いる x = bicgstab (Afcn, b, 1e-6, n, Mfcn)
例5: bicgstab前処理行列M1 とM2 x = bicgstab (A, b, [], 1e-6, n, M1, M2)
例6: bicgstab関数を前処理として用いる x = bicgstab (Afcn, b, 1e-6, n, M1fcn, M2fcn)
例7: bicgstab引数を必要とする関数を入力として
function y = Ap (A, x, z) # compute A^z * x
y = x;
for i = 1:z
y = A * y;
endfor
endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = bicgstab (Apfcn, b, [], [], [], [], [], 2);
例8:bicgstab右前処理子を使用する ことを示す明示的な例 [M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed M = M1 * M2; ## reference solution computed by bicgstab after one iteration [x_ref, fl] = bicgstab (A, b, [], 1, M) ## right preconditioning [y, fl] = bicgstab (A / M, b, [], 1) x = M \ y # compare x and x_ref
参照:
Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM
See also: bicg, cgs, gmres, pcg, qmr, tfqmr.
: x = cgs (A, b, tol, maxit, M1, M2, x0, …)
: x = cgs (A, b, tol, maxit, M, [], x0, …)
: [x, flag, relres, iter, resvec] = cgs (A, b, …)
共役勾配二乗法を使用して、 (Aは正方行列)を解きますA x = b。
入力引数は次のとおりです。
Aは線形システムの行列であり、正方行列でなければなりません。 A はAfcn、のような行列、関数ハンドル、またはインライン関数として渡すことができますAfcn(x) = A * x。 への追加パラメータは、x0Afcnの後に渡されます。
bは右側のベクトルです。これはAと同じ行数の列ベクトルでなければなりません。
tolは相対許容値です。指定されないか [] に設定されている場合、デフォルト値 1e-6 が使用されます。
maxit は外部反復の最大回数です。指定されないか [] に設定されている場合、デフォルト値min (20, numel (b))が使用されます。
M1、M2は前処理行列です。前処理行列は として与えられますM = M1 * M2。M1 とM2 は両方とも、または となるような行列、または関数ハンドルまたはインライン関数として渡すことができます。M1 が空であるか渡されない場合、前処理は適用されません。使用される手法は適切な前処理、 つまりを解いてから です。 gg(x) = M1 \ xg(x) = M2 \ xA*inv(M)*y = bx = inv(M)*y
x0初期推定値。指定されないか [] に設定されている場合、デフォルト値zeros (size (b))が使用されます。
x0に続く引数はパラメータとして扱われ、に渡される関数 ( AまたはPcgs ) のいずれかに適切な方法で渡されます。
出力パラメータは次のとおりです。
vxは計算された近似値です。この方法が収束しない場合は、最小残差で反復されます。
フラグは終了ステータスを示します:
0: 反復は選択された許容範囲内に収束しました
1: 収束前に最大反復回数に達した
2: 前処理行列は特異である
3: アルゴリズムが停滞した
4: ゼロ除算のためアルゴリズムを続行できません
relres は、 として得られた相対残差です 。 (A*x-b) / norm(b)
iter はxが計算される反復です。
resvec は各反復における残差を含むベクトルです。これにより、実行された反復の合計回数を確認することができます。 length(resvec) - 1
三角行列に関する簡単な問題を考えてみましょう
n = 20; A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... sparse (1, 2, 1, 1, n) * n / 2); b = A * ones (n, 1); restart = 5; [M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)' M = M1 * M2; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(x) M2 \ x;
例1:最も簡単な使用法cgs x = cgs (A, b, [], n)
例2: cgs計算する関数 A * x x = cgs (Afcn, b, [], n)
例3: cgs前処理行列M x = cgs (A, b, [], 1e-06, n, M)
例4: cgs関数を前処理として用いる x = cgs (Afcn, b, 1e-6, n, Mfcn)
例5: cgs前処理行列M1 とM2 x = cgs (A, b, [], 1e-6, n, M1, M2)
例6: cgs関数を前処理として用いる x = cgs (Afcn, b, 1e-6, n, M1fcn, M2fcn)
例7: cgs引数を必要とする関数を入力として
function y = Ap (A, x, z) # compute A^z * x
y = x;
for i = 1:z
y = A * y;
endfor
endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = cgs (Apfcn, b, [], [], [], [], [], 2);
例8:cgs右前処理子を使用する ことを示す明示的な例 [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed M = M1 * M2; ## reference solution computed by cgs after one iteration [x_ref, fl] = cgs (A, b, [], 1, M) ## right preconditioning [y, fl] = cgs (A / M, b, [], 1) x = M \ y # compare x and x_ref
References:
Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM
See also: pcg, bicgstab, bicg, gmres, qmr, tfqmr.
: x = gmres (A, b, restart, tol, maxit, M1, M2, x0, …)
: x = gmres (A, b, restart, tol, maxit, M, [], x0, …)
: [x, flag, relres, iter, resvec] = gmres (A, b, …)
A x = b再起動を伴う前処理付き GMRES 反復法 (別名 PGMRES(restart)) を使用して 解きます。
入力引数は次のとおりです。
Aは線形システムの行列であり、正方行列でなければなりません。 A はAfcn、のような行列、関数ハンドル、またはインライン関数として渡すことができますAfcn(x) = A * x。 への追加パラメータは、x0Afcnの後に渡されます。
bは右側のベクトルです。これはAと同じ行数の列ベクトルでなければなりません。
restart は、メソッドが再開するまでの反復回数です。[] または N = numel (b) の場合、再起動は適用されません。
vtol は、前処理された残差誤差に必要な相対許容値です 。 の場合、反復は停止します。tol が省略されるか空の場合、 1e-6 の許容値が使用されます。 inv (M) * (b - a * x)norm (inv (M) * (b - a * x)) ≤ tol * norm (inv (M) * B)maxit は外側の反復の最大回数です。指定されていないか [] に設定されている場合、デフォルト値 が使用されます。restart が空の場合、maxit が反復の最大回数になることに注意してください。restartとmaxitが空でない場合、反復の最大回数は です。restartとmaxit の両方が空の場合、反復の最大回数は に設定されます。 min (10, N / restart)restart * maxitmin (10, N)
M1、M2は前処理子です。前処理子 M はとして与えられますM = M1 * M2。M1とM2 は両方とも、または となるような行列、関数ハンドル、またはインライン関数として渡すことができます。M1 が [] であるか指定されていない場合は、前処理子は適用されません。使用される手法は左前処理子です。つまり、 の代わりに が解かれます。 gg(x) = M1 \ xg(x) = M2 \ xinv(M) * A * x = inv(M) * bA * x = b
x0は初期推定値です。指定されないか [] に設定されている場合、デフォルト値 が使用されます。 zeros (size (b))
x0に続く引数はパラメータとして扱われ、 に渡される関数 ( AまたはMまたは M1またはM2 )のいずれかに適切な方法で渡されますgmres。
出力は次のとおりです。
x は計算された近似値です。この方法が収束しない場合は、最小残差で反復されます。
フラグは終了ステータスを示します:
0 : 反復は指定された許容範囲内に収束しました
1 : 最大反復回数を超えました
2 : 前処理行列は特異である
3:アルゴリズムが停滞した(2つのアルゴリズム間の相対的な差)連続反復回数がeps未満の場合)relres は近似値xの相対的な前処理済み残差の値です。
iter は、 xを計算するために実行された外部反復と内部反復の回数を含むベクトルです。つまり、
iter(1) : 外側の反復回数、つまりメソッドが再起動された回数。( restartが空またはNの場合は1、そうでない場合は1 ≤ iter(1) ≤ maxit )。
iter(2) : 再開前に実行される反復回数。つまり、メソッドは次の場合に再開します 。restartが空または Nの場合、1 ≤ iter(2) ≤ maxitとなります。 iter(2) = restart
より明確に言うと、近似値xは反復 で計算されます 。出力x は最小の事前条件付き残差解に対応するため、この方法で実行される反復の合計数は で与えられます。 (iter(1) - 1) * restart + iter(2)length (resvec) - 1
resvecは、0回目の反復を含む各反復における前処理された相対残差を含むベクトルです 。 norm (A * x0 - b)
三角行列に関する簡単な問題を考えてみましょう
n = 20; A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... sparse (1, 2, 1, 1, n) * n / 2); b = A * ones (n, 1); restart = 5; [M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A) M = M1 * M2; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(x) M2 \ x;
例1:最も簡単な使用法gmres x = gmres (A, b, [], [], n)
例2: gmres計算する関数 A * x x = gmres (Afcn, b, [], [], n)
例3:gmres再起動時の の使用 x = gmres (A, b, restart);
例4: 再起動ありと再起動なしの gmres前処理行列M x = gmres (A, b, [], 1e-06, n, M) x = gmres (A, b, restart, 1e-06, n, M)
例5: gmres関数を前処理として用いる x = gmres (Afcn, b, [], 1e-6, n, Mfcn)
例6: gmres前処理行列M1 とM2 x = gmres (A, b, [], 1e-6, n, M1, M2)
例7: gmres関数を前処理として用いる x = gmres (Afcn, b, 1e-6, n, M1fcn, M2fcn)
例8: gmres引数を必要とする関数を入力として
function y = Ap (A, x, p) # compute A^p * x
y = x;
for i = 1:p
y = A * y;
endfor
endfunction
Apfcn = @(x, p) Ap (A, x, p);
x = gmres (Apfcn, b, [], [], [], [], [], [], 2);
EXAMPLE 9: explicit example to show that gmres uses a left preconditioner [M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed M = M1 * M2; ## reference solution computed by gmres after two iterations [x_ref, fl] = gmres (A, b, [], [], 1, M) ## left preconditioning [x, fl] = gmres (M \ A, M \ b, [], [], 1) x # compare x and x_ref
参照:
Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM
See also: bicg, bicgstab, cgs, pcg, pcr, qmr, tfqmr.
: x = qmr (A, b, rtol, maxit, M1, M2, x0)
: x = qmr (A, b, rtol, maxit, P)
: [x, flag, relres, iter, resvec] = qmr (A, b, …)
A x = b準最小残差反復法(先読みなし)を使用して 解きます。
rtolは相対許容値です。指定されないか [] に設定されている場合、デフォルト値 1e-6 が使用されます。
maxit は外部反復の最大回数です。指定されないか [] に設定されている場合、デフォルト値min (20, numel (b))が使用されます。
vx0初期推定値。指定されないか [] に設定されている場合、デフォルト値zeros (size (b))が使用されます。A は、行列として渡すことも、およびと fなる関数ハンドルまたはインライン関数として渡すこともできます。 f(x, "notransp") = >A*xf(x, "transp") = A'*x
前処理子Pは として与えられますP = M1 * M2。M1 とM2はどちらもまたは かつ またはとなる gような行列または関数ハンドルまたはインライン関数として渡すことができます。 g(x, "notransp") = M1 \ xg(x, "notransp") = M2 \ xg(x, "transp") = M1' \ xg(x, "transp") = M2' \ x
複数の出力パラメータで呼び出された場合
フラグは終了ステータスを示します:
0: 反復は選択された許容範囲内に収束しました
1: 収束前に最大反復回数に達した
3: アルゴリズムが停滞した(値 2 は未使用ですが、互換性のためにスキップされます)。
relres は相対残差の最終値です。
iter は実行される反復回数です。
resvecは、各反復における残差ノルムを含むベクトルです。
参考文献:
R. Freund and N. Nachtigal, QMR: a quasi-minimal residual method for non-Hermitian linear systems, Numerische Mathematik, 1991, 60, pp. 315–339.
R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhour, R. Pozo, C. Romine, and H. van der Vorst, Templates for the solution of linear systems: Building blocks for iterative methods, SIAM, 2nd ed., 1994.
See also: bicg, bicgstab, cgs, gmres, pcg.
: x = tfqmr (A, b, tol, maxit, M1, M2, x0, …)
: x = tfqmr (A, b, tol, maxit, M, [], x0, …)
: [x, flag, relres, iter, resvec] = tfqmr (A, b, …)
A x = bcgs に基づいて、Transpose-Tree qmr 法を使用して 解きます。
入力パラメータは次のとおりです。
Aは線形システムの行列であり、正方行列でなければなりません。 A はAfcn、のような行列、関数ハンドル、またはインライン関数として渡すことができますAfcn(x) = A * x。 への追加パラメータは、x0Afcnの後に渡されます。
vbは右側のベクトルです。Aと同じ行数の列ベクトルでなければなりません。tolは相対許容値です。指定されないか [] に設定されている場合、デフォルト値 1e-6 が使用されます。
maxit は、外側の反復の最大回数です。指定されないか、[] に設定されている場合、デフォルト値がmin (20, numel (b))使用されます。このメソッドは反復回数が奇数か偶数かによって動作が異なるため、互換性を保つために、tfqmr奇数-偶数サイクル全体での反復とみなされます。つまり、反復全体を実行するために、アルゴリズムは 2 つのサブ反復 (奇数反復と偶数反復) を実行します。
M1、M2は前処理子です。前処理子 Mは として与えられますM = M1 * M2。M1とM2 は両方とも、または と なるような行列、または関数ハンドルまたはインライン関数として渡すことができます。使用される手法は右前処理です。つまり、 を解いて から ではなく となり ます。 gg(x) = M1 \ xg(x) = M2 \ xA*inv(M)*y = bx = inv(M)*yA x = b
x0初期推定値。指定されないか [] に設定されている場合、デフォルト値zeros (size (b))が使用されます。
x0に続く引数はパラメータとして扱われ、に渡される関数 ( AまたはMtfqmr ) のいずれかに適切な方法で渡されます。
出力パラメータは次のとおりです。
xは計算された近似値です。この方法が収束しない場合は、最小残差で反復されます。
フラグは終了ステータスを示します:
0: 反復は選択された許容範囲内に収束しました
1: 収束前に最大反復回数に達した
2: 前処理行列は特異である
3: アルゴリズムが停滞した
4: ゼロ除算のためアルゴリズムを続行できません
relres は次のように得られる相対残差です 。 (A*x-b) / norm (b)
iter はxが計算される反復です。
resvec は各反復( を含む)における残差を含むベクトルですnorm (b - A x0)。 を実行すると、実行された反復の合計回数を確認できます。 length (resvec) - 1
三角行列に関する簡単な問題を考えてみましょう
n = 20; A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... sparse (1, 2, 1, 1, n) * n / 2); b = A * ones (n, 1); restart = 5; [M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)' M = M1 * M2; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(x) M2 \ x;
例1:最も簡単な使用法tfqmr x = tfqmr (A, b, [], n)
例2: tfqmr計算する関数 A * x x = tfqmr (Afcn, b, [], n)
例3: tfqmr前処理行列M x = tfqmr (A, b, [], 1e-06, n, M)
例4: tfqmr関数を前処理として用いる x = tfqmr (Afcn, b, 1e-6, n, Mfcn)
例5: tfqmr前処理行列M1 とM2 x = tfqmr (A, b, [], 1e-6, n, M1, M2)
例6: tfmqr関数を前処理として用いる x = tfqmr (Afcn, b, 1e-6, n, M1fcn, M2fcn)
例7: tfqmr引数を必要とする関数を入力として
function y = Ap (A, x, z) # compute A^z * x
y = x;
for i = 1:z
y = A * y;
endfor
endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = tfqmr (Apfcn, b, [], [], [], [], [], 2);
EXAMPLE 8: explicit example to show that tfqmr uses a right preconditioner [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed M = M1 * M2; ## reference solution computed by tfqmr after one iteration [x_ref, fl] = tfqmr (A, b, [], 1, M) ## right preconditioning [y, fl] = tfqmr (A / M, b, [], 1) x = M \ y # compare x and x_ref
参照:Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM
See also: bicg, bicgstab, cgs, gmres, pcg, qmr, pcr.