ek_struct.m

function [S, ST] = ek_struct(A, cholesky)
%BUILD_RK_STRUCT Build a struct for building rational Krylov subspaces.
%

if issparse(A)
    if exist('cholesky', 'var') && cholesky
        [R, g, pA] = chol(A);
        
        if (g ~= 0)
            % warning('Matrix is not posdef, resorting to LU');
            [S, ST] = ek_struct(A, false);
            return;
        end
        
        S = struct(...
            'solve', @(nu, mu, x) sparse_solve(nu, mu, R', R, pA', pA, x), ...
            'multiply', @(rho, eta, x) rho * A * x - eta * x, ...
            'isreal', isreal(A), ...
            'nrm', normest(A, 1e-2));
        ST = S;
        
    else
        [LA, UA, pA, qA] = lu(A);
        
        S = struct(...
            'solve', @(nu, mu, x) sparse_solve(nu, mu, LA, UA, pA, qA, x), ...
            'multiply', @(rho, eta, x) rho * A * x - eta * x, ...
            'isreal', isreal(A), ...
            'nrm', normest(A, 1e-2));
        
        if nargout >= 2
            ST = struct(...
                'solve', @(nu, mu, x) sparse_solve(nu, mu, UA', LA', pA', qA', x), ...
                'multiply', @(rho, eta, x) rho * A' * x - eta * x, ...
                'isreal', S.isreal, ...
                'nrm', S.nrm);
        end
    end
    
elseif isa(A, 'hm')
    [LA,UA] = lu(A);
    
    S = struct(...
        ... % 'solve', @(nu, mu, x) (nu * A - mu * eye(size(A), 'like', A)) \ x, ...
        'solve', @(nu, mu, x) lu_solve(nu, mu, LA, UA, x),...
        'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
        'isreal', isreal(A), ...
        'nrm', norm(A));
    
    if nargout >= 2
        ST = struct(...
            ...% 'solve', @(nu, mu, x) (nu * A' - mu * eye(size(A), 'like', A)) \ x, ...
            'solve', @(nu, mu, x) lu_solve(nu, mu, UA', LA', x), ...
            'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
            'isreal', S.isreal, ...
            'nrm', S.nrm);
    end
    
elseif isa(A, 'hss')
    ULV = ulv(A);
    
    S = struct(...
        ... % 'solve', @(nu, mu, x) (nu * A - mu * eye(size(A), 'like', A)) \ x, ...
        'solve', @(nu, mu, x) ulv_solve_wrapper(nu, mu, ULV, x),...
        'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
        'isreal', isreal(A), ...
        'nrm', norm(A));
    
    if nargout >= 2
        ULVT = ulv(A');
        ST = struct(...
            ...% 'solve', @(nu, mu, x) (nu * A' - mu * eye(size(A), 'like', A)) \ x, ...
            'solve', @(nu, mu, x) ulv_solve_wrapper(nu, mu, ULVT, x), ...
            'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
            'isreal', S.isreal, ...
            'nrm', S.nrm);
    end
end

end

function y = sparse_solve(nu, mu, L, U, p, q, x)

if nu > mu
    y = nu\(q*(U\(L\(p*x))));
else
    y = -mu \ x;
end

end