% Experiment 1 - first part

clear all;
close all;

n_vec = [600, 800 1000];
data_plot = zeros(29, 6);
t = 1;

for i = 1:length(n_vec)
    n = n_vec(i);
    e = ones(n, 1);
    A = spdiags([e -2 * e e], -1:1, n, n);
    I = speye(n);
    A = kron(A, I) + kron(I, A);

    for q = [1 3 5]
        rng('default');
        B = randn(n^2, q);
        B = B / norm(B, 'fro');

        opts.maxit_m = 100;
        opts.maxit_j = 500;
        opts.tol_outer = 1e-8;
        opts.tol_inner = 1e-10;
        opts.type = 'Galerkin';

        % LR-ADI-EKSM routine
        tt = tic;
        [Z, inner_dimension, pvec, lyap_residual, time] = LR_ADI_EKSM(A, B, opts);
        time_LR_ADI_EKSM = toc(tt);

        semilogy(lyap_residual,'+-')
        pause
        
        % compute the real residual norm
        fact = [A * Z, Z, B];
        [~, R] = qr(fact, 0);
        Theta = zeros(size(fact, 2));
        Theta(1:size(Z, 2), size(Z, 2) + 1:2 * size(Z, 2)) = eye(size(Z, 2));
        Theta(size(Z, 2) + 1:2 * size(Z, 2), 1:size(Z, 2)) = eye(size(Z, 2));
        Theta(2 * size(Z, 2) + 1:end, 2 * size(Z, 2) + 1:end) = eye(size(B, 2));
        res_real = norm(R * Theta * R', 'fro') / norm(B' * B, 'fro');
        fprintf(['LR-ADI-EKSM, n: %6d, q: %3d\n ADI its: %3d\n' ...
                 'Times\n Basis: %10.5e, Prj. Pr.: %10.5e, Shifts: ' ...
                 '%10.5e, Etc: %10.5e, Tot: %10.5e\n'], n^2, q, ...
                length(lyap_residual), time.basis, time.projected_problem, ...
                time.shifts, time.etc, time_LR_ADI_EKSM);

        if i == 1
            data_plot(1:length(lyap_residual), t) = 1:length(lyap_residual);
            data_plot(1:length(lyap_residual), t + 1) = inner_dimension;
            data_plot(1:length(lyap_residual), t + 2) = lyap_residual;
            
            t = t + 2;
        end
    end
end
writematrix(data_plot, ...
            'results_Experiment1_firstpart.txt', ...
            'Delimiter', 'tab');
