% Expweriment 2

clear all;
warning off;
addpath(genpath('./'));
n_vec = [50, 80, 100];
% n_vec=50
data_plot = zeros(80, 4);
r = 1;
% epsilon: viscosity parameter. For epsilon << 1 the problem becomes harder
% as we have dominant convection
epsilon_vec = [.05, .005];

for i = 1:length(epsilon_vec)
    epsilon = epsilon_vec(i);

    for j = 1:length(n_vec)
        n = n_vec(j);

        h = 1 / (n - 1);

        % Negative Laplacian
        T = -epsilon / h^2 * ...
            spdiags([-ones(n, 1), 2 * ones(n, 1), -ones(n, 1)], -1:1, n, n);
        % negative first derivative
        N = -1 / (2 * h) * ...
            spdiags([-ones(n, 1), zeros(n, 1), ones(n, 1)], -1:1, n, n);

        x = linspace(0, 1, n);

        phi1 = 1 - x.^2;
        psi1 = x;
        upsilon1 = x;
        upsilon3 = exp(x);

        PHI1 = spdiags(phi1', 0, n, n);
        PSI1 = spdiags(psi1', 0, n, n);
        UPSILON1 = spdiags(upsilon1', 0, n, n);
        UPSILON3 = spdiags(upsilon3', 0, n, n);

        B1 = PHI1 * N;
        B3 = UPSILON3 * N';

        I = speye(n);
        % coefficient matrix
        A = kron(kron(I, T), I) + ...
            kron(T + B3, kron(I, I)) + ...
            kron(kron(I, I), T) + ...
            kron(UPSILON1, kron(PSI1, B1));
        % preconditioner
        ups1 = mean(upsilon1);
        psi1 = mean(psi1);
        opts.P1 = full(T + B3);
        opts.P2 = full(kron(I, T) + kron(T, I) + kron(ups1 * PSI1, B1));

        % opts.P=1;

        rng('default');
        B = randn(n^3, r);

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

        tt = tic;
        [Z, inner_dimension, p_vec, lyap_residual, time, inner_it_GMRES] = ...
            LR_ADI_EKSM(A, B, opts);
        time_LR_ADI_EKSM = toc(tt);

        fact = [A * Z, Z, B];
        [~, R] = qr(fact, 0);
        Theta = zeros(size(fact, 2));
        one = 1:size(Z, 2);
        two = size(Z, 2) + 1:2 * size(Z, 2);
        three = 2 * size(Z, 2) + 1:size(fact, 2);
        Theta(one, two) = eye(size(Z, 2));
        Theta(two, one) = eye(size(Z, 2));
        Theta(three, three) = eye(size(B, 2));
        res_real = norm(R * Theta * R', 'fro') / norm(B' * B, 'fro');

        complex_shifts = length(find(imag(p_vec) ~= 0));

        fprintf(['Problem dimension: %6d, epsilon=%10.5e, #(pure ' ...
                 'complex shifts): %5d\n'], ...
                n^3, ...
                epsilon, ...
                complex_shifts);
        fprintf(['LR-ADI-EKSM\n ADI its: %3d, Time: %10.5e, Rel. ' ...
                 'Res: %10.5e\n'], ...
                length(lyap_residual), ...
                time_LR_ADI_EKSM - time.shifts, ...
                res_real);

        % standard LR-ADI routine
        % use same shifts
        opts.p = p_vec;
        opts.prec_shifted = 1;
        tt = tic;
        [Z_lradi, lyap_residual2] = LR_ADI(A, B, opts);
        time_LR_ADI = toc(tt);

        fact = [A * Z_lradi, Z_lradi, B];
        [~, R] = qr(fact, 0);
        Theta = zeros(size(fact, 2));
        one = 1:size(Z_lradi, 2);
        two = size(Z_lradi, 2) + 1:2 * size(Z_lradi, 2);
        three = 2 * size(Z_lradi, 2) + 1:size(fact, 2);
        Theta(one, two) = eye(size(Z_lradi, 2));
        Theta(two, one) = eye(size(Z_lradi, 2));
        Theta(three, three) = eye(size(B, 2));
        res_real = norm(R * Theta * R', 'fro') / norm(B' * B, 'fro');

        fprintf(['LR-ADI\n ADI its: %3d, Time: %10.5e, Rel. Res: ' ...
                 '%10.5e\n'], ...
                length(lyap_residual2), ...
                time_LR_ADI, ...
                res_real);

        %%% ADI MMESS
        mess_path();
        eqn.A_ = A;
        eqn.B = B;
        eqn.P1 = opts.P1;
        eqn.P2 = opts.P2;

        opts.prec_shifted = 1;
        eqn.type = 'N';
        eqn.haveE = 0;
        eqn.haveUV = 0;
        opts.adi.res_tol = 1e-8;
        opts.shifts.p = p_vec;
        % set operation
        oper = operatormanager('default');
        % change oper.sol_ApE to allow for an iterative solution of the
        % shifted linear systems
        oper.sol_ApE = @iterative_sol;

        tt = tic;
        [out, eqn, opts, oper] = mess_lradi(eqn, opts, oper);
        time_LR_ADI_MMESS = toc(tt);
        Z_mmess = out.Z;

        fact = [A * Z_mmess, Z_mmess, B];
        [~, R] = qr(fact, 0);
        Theta = zeros(size(fact, 2));
        one = 1:size(Z_mmess, 2);
        two = size(Z_mmess, 2) + 1:2 * size(Z_mmess, 2);
        three = 2 * size(Z_mmess, 2) + 1:size(fact, 2);
        Theta(one, two) = eye(size(Z, 2));
        Theta(two, one) = eye(size(Z, 2));
        Theta(three, three) = eye(size(B, 2));
        res_real = norm(R * Theta * R', 'fro') / norm(B' * B, 'fro');

        fprintf(['LR-ADI_MMESS\n ADI its: %3d, '...
                 'Time: %10.5e, Rel. Res: %10.5e\n'], ...
                out.niter, ...
                time_LR_ADI_MMESS, ...
                res_real);

        if n <= 50
            myrange = 1:length(inner_it_GMRES);
            data_plot(myrange, (i - 1) * 2 + 1) = myrange;
            data_plot(myrange, (i - 1) * 2 + 2) = inner_it_GMRES;
        end

        clear p_vec;
        clear opts;
        clear eqn;

    end
end
writematrix(data_plot, 'results_Experiment2.txt', 'Delimiter', 'tab');

