% 連立一次方程式に対する CG 法（MATLABルーチン使用版）サンプル
clear all;
clc;
close all;

% 諸パラメータ（ここをいろいろ変更してみるとよい）
n = 400;        % 行列のサイズ（※下記で若干変更が入る場合があるので注意）
tol = 1e-12;

% 問題（係数行列）の定義（いろいろ試してみるとよい）
prob = 1;           % 問題を切り替え
switch prob
    case 1
        m = round(sqrt(n));
        A = gallery('poisson',m);    % 2次元Laplacianの離散化行列（→Poisson方程式）
        n = m * m;
    case 2
        A = sparse(gallery('toeppd', n));    % Toeplitz行列（密行列）
    case 3
        m = round(sqrt(n));
        A = sparse(gallery('wathen', m, m));    % 有限要素法の行列らしい
        [n m] = size(A);
    case 4
        A = sparse(gallery('ris', n) + 2 * eye(n) );  % 対称Hankel行列 + 2I　（密行列）
    case 5
        % SuitSparseMatrixCollection のデータを使う場合はこのようにする
        % あらかじめ mat ファイルをダウンロードしてソースと同じところに置いておくこと
        load('1138_bus.mat');   % 例：1138_bus（疎行列）
        A = Problem.A;
        [n m] = size(A);
end


% 前処理行列の定義
%opts.type='ict';    % ICのタイプ
opts.type='nofill';    % ICのタイプ
M = ichol(A);       % IC(0)

% 問題（右辺ベクトル）の定義
x = rand(n,1);      % サンプルの解
b = A*x;            % 問題を作る Ax = b

% 前処理なしのCG法
[xx1, fl1, rr1, it1, rv1] = pcg(A, b, tol, 2*n);

% 前処理付きのCG法
[xx2, fl2, rr2, it2, rv2] = pcg(A, b, tol, 2*n, M, M');

% 図を描く（このブロックは MATLAB ドキュメントからのほぼ引用です）
figure;
semilogy(rv1/norm(b), 'b.');
hold on;
semilogy(rv2/norm(b), 'r.');
legend('No Preconditioner','IC(0)');
title(['Dimension n = ', num2str(n)]);
xlabel('iteration number');
ylabel('relative residual');
hold off;

% 条件数の推定値を計算
[condest(A), condest(M'\(M\A))]