% 固有値問題に対する QR 法（原型版）サンプル
% 　下記アルゴリズムは推奨される実装ではなくあくまで動作サンプルです
%   現実の実装から乖離している（不足している）点：
%   　- 反復が高速化されておらず O(n^3) かかる（ので Hessenberg 化は本来の効力を発揮してない）
%     - 収束高速化（シフト導入，減次による問題縮減）が組み込まれていない（講義未説明）
%     - 複素固有値ペアの扱いがなされていない（講義未説明）
%   したがって反復はかなり遅い（計算量的にも反復回数的にも）です
clear all;
clc;
close all;

% MATLAB の出力（表示）桁数を倍精度にするおまじない
%format long;

% 諸パラメータ（ここをいろいろ変更してみるとよい）
tol = 1e-12;        % 停止則に使う相対誤差
max_itr = 10000;    % 強制打ち切りの反復回数（収束状況で加減して様子を見てください）
snapshot = 0;       % 1 = 定期的に行列を出力する，それ以外 = 行列を一切表示しない
every = 100;        % 上が 1 の場合，every 反復ごとに行列を出力する

% 問題（係数行列）の定義（いろいろ試してみるとよい）
prob = 1;           % 問題を切り替え
switch prob
    case 1      % 自分で問題を設定するやり方
        % 指定固有値を持つランダム行列（ベクトル v に任意の固有値を並べる）
        v = 1:10;                   % 行列の固有値をこれで指定
        %v = [v -1 -2 -3];           % 難しくするには絶対値の等しい固有値を追加する
        A = diag(v);                % 対角行列
        [n m] = size(A);            % A の大きさを取得
        S = rand(n);                % ランダム行列（正則だと期待）
        A = S * A * inv(S);         % それで A を相似変換（逆行列は使うなと言いましたが）
    case 2      % とても難しい例
        % 次の行列は MATLAB のテスト行列に入っている 5x5 の実行列
        % - eig(A) してみると複素固有値ペアが2組ある（のでそいつらは 2x2 ブロックに行く）
        % - しかも abs(eig(A)) してみるとすべての固有値がほとんど大きさが同じ
        % ので，この問題は，5x5の割には難しい
        % こういう問題を解くにはさらに工夫が必要（講義では説明していない各種技法が必要）
        A = gallery(5);
end

% Hessenberg 化
H = hess(A);        % MATLAB の Hessenberg 化ルーチンにお願い
Hnorm = norm(H);    % H のノルムを計算しておく

% ナイーブな QR 法の反復
for i = 1:max_itr
    [Q R] = qr(H);     % MATLAB の QR 分解ルーチンにお願い
    H = R*Q;
    if subd(H)/Hnorm < tol   % 副対角が十分小さくなったらループを抜ける
        break;
    end
    if snapshot == 1     % スナップショットを表示する場合の手続き
        if mod(i, every) == 0
            i           % 反復数表示
            H           % その時点の行列 H を表示
            pause;      % 何かキーを押すまで一次停止
        end
    end
end
i                   % 反復数

% お答を比較
e1 = sort(diag(H), 'descend');       % 上のナイーブ QR 法のお答え
e2 = sort(eig(A), 'descend');   % MATLAB の固有値ルーチンのお答え

[e1 e2]             % 並べて眺めてみる
norm(e1-e2)         % ノルムを計ってみたりとか


% 以降，途中で使われるサブルーチン
function res = subd(A)
% 入力行列の下副対角の絶対値最大の値を返す
    [n m] = size(A);
    res = 0.0;
    for i=1:n-1
        res = max( abs(A(i+1,i)), res );    % もっと賢い方法はないものか
    end
end