% Experiment 3
clear all;
load MAX_Gestell_matrices.mat;
E_cell = E;
B_cell = B;
C_cell = C;
addpath(genpath('./'));

for k = [4 3 1]
    EE = EGam{k};
    A = tA{k};
    for i = 1:length(EE)
        A = A - 10 * (EE{i});
    end
    E = E_cell{k};
    % the right-hand side for k=3 is rank deficient:
    % it has rank 5 even though it has 6 columns
    B = B_cell{k};
    C = C_cell{k};
    if k == 3
        [Q, R] = qr(full(B), 0);
        [U, S, V] = svd(full(R), 0);
        B = Q * (U(:, 1:5) * S(1:5, 1:5));
    end
    fprintf('n: %6d, rank(rhs): %6d\n', size(A, 1), size(B, 2));

    n = size(A, 1);
    Esqrtinv = spdiags(sqrt(1 ./ diag(E)), 0:0, n, n);
    eqn.A_ = Esqrtinv * A * Esqrtinv;
    eqn.B = Esqrtinv * B;

    eqn.C = C;
    eqn.E_ = speye(size(A));
    eqn.type = 'N';
    eqn.haveE = 1;
    eqn.haveUV = 0;

    opts.maxit_m = 100;
    opts.maxit_j = 100;
    opts.tol_inner = 1e-10;
    opts.tol_outer = 1e-8;

    opts.type = 'Galerkin';
    opts_eigs.tols = 1e-2;
    a = eigs(eqn.A_, 1, 'sm', opts_eigs);
    b = eigs(eqn.A_, 1, 'lm', opts_eigs);
    tol = opts.tol_outer * norm(eqn.B' * eqn.B, 'fro') / 1e1;
    opts.p = myadipars(a, b, 0, tol);

    tt = tic;
    [Z, inner_dimension, p_vec, lyap_residual, time] = ...
        LR_ADI_EKSM(eqn.A_, full(eqn.B), opts);
    time_LR_ADI_EKSM = toc(tt);
    fprintf('LR-ADI-EKSM\n ADI its: %3d, Time: %10.5e, dim(EK): %5d\n', ...
            length(lyap_residual), ...
            time_LR_ADI_EKSM - time.shifts, ...
            inner_dimension(end));

    mess_path(1);

    % set operation
    oper = operatormanager('state_space_transformed_default');

    eqn.type = 'N';
    opts.KSM.maxiter = 500;
    opts.KSM.restol = 1e-8;
    opts.KSM.space = 'EK';
    opts.KSM.type_eqn = 'LE';
    opts.KSM.LDL_T = 0;
    opts.KSM.symmetric = 1;
    opts.KSM.tol_trunc = 1e-16;
    opts.KSM.explicit_proj = 0;
    opts.KSM.rctol = 0;
    opts.KSM.comp_res = 1;

    opts.KSM.info = 0;

    fprintf('\nComputing solution factorization.\n');

    tt = tic;
    [out, eqn, opts, oper] = mess_KSM(eqn, opts, oper);
    time_KPIK = toc(tt);

    fprintf('KPIK\n ADI its: %3d, Time: %10.5e, dim(EK): %5d\n', ...
            out.it, time_KPIK, 2 * (out.it + 1) * size(B, 2));

    Z_KPIK = Esqrtinv * out.Z;

    %%% ADI MMESS

    opts.adi.res_tol = 1e-8;
    opts.adi.maxiter = 500;
    opts.adi.compute_sol_fac = 1;
    opts.adi.info = 0;
    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 = Esqrtinv * out.Z;

    fprintf('LR-ADI_MMESS\n ADI its: %3d, Time: %10.5e\n', ...
            out.niter, time_mmess);

    clear p_vec;
    clear opts;
    clear eqn;

end
