Keyword: スパース, 固有値問題
概要
本サンプルはスパース固有値問題を解くサンプルプログラムです。 本サンプルは以下に示されるスパース固有値問題を解いて固有値を出力します。本サンプルプログラムでは、関数 f12ac のほか、f12aa、f12ab、f12ad や f12ae を呼び出しています。
※本サンプルはnAG Toolbox for MATLAB®が提供する関数 f12ac() のExampleコードです。実行にはMATLAB®本体（他社製品）とnAG Toolbox for MATLAB®が必要です。
本サンプル及び関数の詳細情報は f12ac のマニュアルページをご参照ください。
入力データ
n = int32(100);
nev = int32(4);
ncv = int32(20);
h = 1/(double(n)+1);
rho = 10;
md = repmat(4*h, double(n), 1);
me = repmat(h, double(n-1), 1);
irevcm = int32(0);
resid = zeros(n,1);
v = zeros(n, ncv);
x = zeros(n, 1);
mx = zeros(n);
dd = 2/h;
dl = -1/h - rho/2;
du = -1/h + rho/2;
y = zeros(n,1);
[icomm, comm, ifail] = f12aa(n, nev, ncv);
[icomm, comm, ifail] = f12ad('REGULAR INVERSE', icomm, comm);
[icomm, comm, ifail] = f12ad('GENERALIZED', icomm, comm);
% Construct m and factorise
[md, me, info] = f07jd(md, me);
while (irevcm ~= 5)
[irevcm, resid, v, x, mx, nshift, comm, icomm, ifail] = ...
f12ab(irevcm, resid, v, x, mx, comm, icomm);
if (irevcm == -1 || irevcm == 1)
y(1) = dd*x(1) + du*x(2);
for i = 2:n-1
y(i) = dl*x(i-1) + dd*x(i) + du*x(i+1);
end
y(n) = dl*x(n-1) + dd*x(n);
[x, info] = f07je(md, me, y);
elseif (irevcm == 2)
y(1) = 4*x(1) + x(2);
for i=2:n-1
y(i) = x(i-1) + 4*x(i) + x(i+1);
end
y(n) = x(n-1) + 4*x(n);
x = h*y;
elseif (irevcm == 4)
[niter, nconv, ritzr, ritzi, rzest] = f12ae(icomm, comm);
if (niter == 1)
fprintf('\n');
end
fprintf('Iteration %2d No. converged = %d Norm of estimates = %16.8e\n', niter, nconv, norm(rzest));
end
end
[nconv, dr, di, z, v, comm, icomm, ifail] = f12ac(0, 0, resid, v, comm, icomm)
- １〜３行目は f12aa の入力パラメータを指定しています。
- ５〜６行目は計算に使用するための h と ρ（ロー）の値を指定しています。
- ７〜８行目は f07je の入力パラメータを指定しています。
- １０〜１４行目は f12ab の入力パラメータを指定しています。
- １６〜１９行目は計算に使用する dd、dl、du、yの値を指定しています。
- ２１行目は初期化のため f12aa を呼び出しています。
- ２２行目はオプション設定のため f12ad を呼び出しています。
- ２３行目はオプション設定のため f12ad を呼び出しています。
- ２６行目は実対称正定値三重対角連立一次方程式の計算のため f07jd を呼び出しています。
- ２８〜５２行目は反復して f12ab を呼び出し固有値を求めています。
- 最後に本関数 f12ac を呼び出す構文を指定しています。
出力結果
Iteration 1 No. converged = 0 Norm of estimates = 5.56325675e+03
Iteration 2 No. converged = 0 Norm of estimates = 5.44836588e+03
Iteration 3 No. converged = 0 Norm of estimates = 5.30320774e+03
Iteration 4 No. converged = 0 Norm of estimates = 6.24234186e+03
Iteration 5 No. converged = 0 Norm of estimates = 7.15674705e+03
Iteration 6 No. converged = 0 Norm of estimates = 5.45460864e+03
Iteration 7 No. converged = 0 Norm of estimates = 6.43147571e+03
Iteration 8 No. converged = 0 Norm of estimates = 5.11241161e+03
Iteration 9 No. converged = 0 Norm of estimates = 7.19327824e+03
Iteration 10 No. converged = 1 Norm of estimates = 5.77945489e+03
Iteration 11 No. converged = 2 Norm of estimates = 4.73125738e+03
Iteration 12 No. converged = 3 Norm of estimates = 5.00078500e+03
nconv =
4
dr =
1.0e+04 *
2.0383
2.0339
2.0265
2.0163
di =
0
0
0
0
z =
array elided
v =
array elided
comm =
array elided
icomm =
array elided
ifail =
0
- １〜１２行目は各反復の回数、収束回数とノルム推定を示しています。
- nconv は収束した固有値の数を示しています。
- dr は収束した近似固有値の実部を示しています。
- di は収束した近似固有値の虚部を示しています。
- z は dr と di に対応する固有ベクトルを示します。ここでは出力が省略されています。
- v は近似シュール（Schur）ベクトルを示します。ここでは出力が省略されています。
- comm は現在の解のデータを示します。ここでは出力が省略されています。
- icomm は現在の解のデータを示します。ここでは出力が省略されています。
- ifail は関数がエラーを検知しなければ"0"を出力します。