関連情報

スパース固有値問題

NAG Toolbox for MATLAB®用のサンプルコード

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)

  • 1〜3行目は f12aa の入力パラメータを指定しています。
  • 5〜6行目は計算に使用するための h と ρ(ロー)の値を指定しています。
  • 7〜8行目は f07je の入力パラメータを指定しています。 
  • 10〜14行目は f12ab の入力パラメータを指定しています。
  • 16〜19行目は計算に使用する dd、dl、du、yの値を指定しています。
  • 21行目は初期化のため f12aa を呼び出しています。
  • 22行目はオプション設定のため f12ad を呼び出しています。
  • 23行目はオプション設定のため f12ad を呼び出しています。
  • 26行目は実対称正定値三重対角連立一次方程式の計算のため f07jd を呼び出しています。
  • 28〜52行目は反復して 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

  • 1〜12行目は各反復の回数、収束回数とノルム推定を示しています。
  • nconv は収束した固有値の数を示しています。
  • dr は収束した近似固有値の実部を示しています。
  • di は収束した近似固有値の虚部を示しています。
  • z は dr と di に対応する固有ベクトルを示します。ここでは出力が省略されています。
  • v は近似シュール(Schur)ベクトルを示します。ここでは出力が省略されています。
  • comm は現在の解のデータを示します。ここでは出力が省略されています。
  • icomm は現在の解のデータを示します。ここでは出力が省略されています。
  • ifail は関数がエラーを検知しなければ"0"を出力します。

Results matter. Trust NAG.

Privacy Policy | Trademarks