% Experiment 4

clear all;
addpath(genpath('./'));
mess_path();

for k = 3:4

    eqn = mess_get_linear_rail(k);

    A = eqn.A_;
    E = eqn.E_;
    B = eqn.B;
    C = eqn.C;

    opts.maxit_m = 500;
    opts.maxit_j = 100;
    opts.tol_outer = 1e-8;
    opts.tol_inner = 1e-10;
    pp = symamd(E);
    E_EKSM = E(pp, pp);
    A_EKSM = A(pp, pp);
    C_EKSM = C(:, pp);

    opts.L = chol(E_EKSM, 'lower');
    opts.type = 'Galerkin';
    opts_eigs.tols = 1e-2;

    tt = tic;
    [Z, inner_dimension, p_vec, lyap_residual, time] = ...
        LR_ADI_EKSM(A_EKSM', full(C_EKSM'), opts);
    time_LR_ADI_EKSM = toc(tt);
    Z_LR_ADI_EKSM(pp, :) = Z;

    fprintf('Problem dimension: %6d\n', size(A, 1));
    fprintf(['LR-ADI-EKSM\n ADI its: %3d, Time: %10.5e, dim(EK): ' ...
             '%5d, rank(X): %5d\n'], ...
            length(lyap_residual), time_LR_ADI_EKSM - time.shifts, ...
            inner_dimension(end), size(Z_LR_ADI_EKSM, 2));

    %%% ADI MMESS

    eqn.haveE = 1;
    eqn.type = 'T';
    opts.adi.res_tol = 1e-8;
    opts.adi.maxiter = 500;
    opts.adi.compute_sol_fac = 1;
    opts.adi.info = 1;
    opts.shifts.p = p_vec;
    opts.shifts.method = 'projection';
    opts.shifts.info = 0;
    % set operation
    clear oper;
    oper = operatormanager('default');

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

    fprintf(['LR-ADI_MMESS\n ADI its: %3d, Time: %10.5e, rank(X): ' ...
             '%5d\n'], ...
            out.niter, time_mmess, size(Z_mmess, 2));

    clear Z_LR_ADI_EKSM;

end
