function [funAb,xi,err,V,Z,Y,beta] = RKSM_TLBT_ind1(A,E,B,st,sfreq, tol,maxit,sh,npts,T,ssolve,modG,opts)
% Rational Krylov method for computing low-rank solution factors of the
% solution of the time-limited Gramian equation
%   SXM'+MXS-BtBt'+FF'=0, where S,M,F are index-0 representation
%    of index-1(A,E,B) DAE 
%       and Bt = exp(M\S*T)*B
%
%  CALLING SEQUENCE:
%     [funAb,xi,err,V,Z,Y,beta] = RKSM_TLBT_ind1(A,E,v,st,sfreq, tol,maxit,sh,exact,npts,T,ssolve,modG,opts)
%
% INPUT:
%    E,A    -     n-x-n matrices
%    B      -     n-x-m matrix (m << n);
%    st     - number of dynamic states   
%    sfreq  -    reduced expm(..)x and GCALE is computed/solved every sfreq iteration
%    tol    -     stopping tolerance(s) w.r.t. 
%                   tol(1): y=f(A)B, tol(2) GCALE residual (default tol(1)), 
%                   and tol(3) rank truncation in the end (default 0)
%    maxit  -     maximal number iterations

%    sh     -    shifts (poles) for RKSM (default emptry-> adaptive shifts)
%               can also be a string for the adaptive shifts:
%               'conv' -- choose from convex hull of Ritz values
%               'convR' -- real version of 'conv'
%               'real' -- uses predetermined real interval (outdated)
%    exact  -    exact y for testing        
%    npts   -     numbers of test points for adaptive shift selection
%    T    -     finite time value(s)
%    ssolve  - type of small f(A) handling, currently not used
%    modG    - set to1 for modified time-limited Gramians (rhs is transformed [B,y]) 
% opts    = structure containing info on mass matrix handling 
%           opts.ML,opts.MU --Cholesky or LU factors of E=ML*MU, if E=E^T,
%           MU=ML'
%           opts.facM - compute factors of E if empty ML,MU, 
%               facM=1: compute Chol. / LU factors of E 
%               facM=2: compute Chol. / LU factors of E with pivoting, DEFAULT  
%               facM=0: or directly work with E (facM=0)

% OUTPUT:
%    funAb      approx of exp(A*T)B
%    xi         generated / used shifts   
%    err       vector containing the relative changes
%               ||ynew-yold||/||ynew||, uses exact if given
%    V         orthonormal basis matrix of the rational Krylov subspace
%    Z,Y       low-rank solution factors of GCALE solution ZYZ'~X
%    beta      orthogonalization factors of B 
%
%-----------------------------------------------------------------------
% REFERENCE:
% P. Kuerschner, Balanced truncation model order reduction in 
% limited time intervals for large systems, 
% Advances in Computational Mathematics, 44, pp 1821–1844, 2018;
% https://doi.org/10.1007/s10444-018-9608-6
% 
% Copyright (C) 2015-22 Patrick Kürschner
% 
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1. Redistributions of source code must retain the above copyright notice,
%    this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% INPUT CHECKING
%
if nargin < 13 || isempty(opts),
    opts.facM=0; opts.ML=[]; opts.MU=[];
    opts.saveLU_A=2;
end

if nargin < 10 || isempty(npts), npts = 1601; end
if nargin < 8 || isempty(sh),
    sh='conv'; adaptive = 1;%xi_cand = 1i*logspace(log10(f(1)),log10(f(2)),npts);
end

if size(tol,2)>2 % stopping tolerances
    tolP=tol(2);
    tolsol=tol(3);
    tol=tol(1);
else tolP=tol; 
    tolsol=0;
end
if nargin < 12 || isempty(modG), modG=0; end

if nargin<11 || isempty(ssolve) % currently not used
    ssolve='expm';
end

if (norm(A-A',1)<1e-14), symA=1; else symA=0;end

if isnumeric(sh)
    xi=sh;adaptive = 0;
else
    switch(sh)
        case 'conv'
            adaptive = 1; xi=[];
        case 'convR'
            adaptive = 1; xi=[];
            xi_cand = logspace(-10,10,npts);
    end;
end
[N,s] = size(B);
timings=zeros(5,maxit+1);
timings(1,1)=0;
% handling of mass matrix
pM=1:st;
if isempty(E)
    ML=speye(st); MU=ML; E=spdiags([ones(st,1);zeros(N-st,1)],0,N,N);
else
    if (norm(E-E',1)<1e-14), symM=1; else symM=0;end
    if ~isempty(opts.ML)
        ML=opts.ML;
        if ~isempty(opts.MU), MU=opts.MU;
        else
            if symM, MU=ML'; else MU=speye(N); end
        end
    else
        switch opts.facM
            case 1 % use LU / Cholesky factorization
                ls=tic;
                if symM
                    ML=chol(E(1:st,1:st),'lower'); MU=ML';
                else
                    [ML,MU]=lu(E(1:st,1:st));
                end
                timings(1,1)=toc(ls);
            case 2 % use LU / Cholesky factorization with pivoting
                ls=tic;
                if symM
                    [ML,~,pM]=chol(E(1:st,1:st),'lower','vector'); MU=ML';
                    qM(pM)=1:st;
                else
                    [ML,MU,pM,qM]=lu(E(1:st,1:st),'vector');
                end
                timings(1,1)=toc(ls);
            otherwise %no factor.
                ML=E(1:st,1:st);
                MU=speye(st);
        end
    end
end
switch opts.saveLU_A
    case 1  % use LU factorization
        if symA
            [LA]=chol(-A(st+1:N,st+1:N),'lower'); UA=-LA';
        else
            [LA,UA]=lu(A(st+1:N,st+1:N));
        end
    case 2 % use LU factorization with pivoting
        if symA
            [LA,~,p]=chol(-A(st+1:N,st+1:N),'lower','vector'); UA=-LA';
            q=p;
        else
            [LA,UA,p,q]=lu(A(st+1:N,st+1:N),'vector');
        end
end
if modG
    flip=blkdiag(eye(s),-eye(s));
    modG_rhs=0;
end
if length(T)>1, if min(T)==0, T=max(T); end, end
exact=[];
if ~isempty(exact), 
    switch opts.facM,
        case 2
            exact=full(ML\exact(pM,:));
        otherwise            
            exact=full(ML\exact);
    end
    Nexact=norm(exact); 
end
Is=speye(s);nrmrestot=[];
ls=tic; B=full(B);
switch opts.saveLU_A
    case 1  % use LU factorization
        w1=UA\(LA\B(st+1:N,:));
    case 2 % use LU factorization with pivoting
        rhs=B(st+1:N,:);
        w1(q,:)=UA\(LA\rhs(p,:));
    otherwise %no factor.
        w1=A(st+1:N,st+1:N)\B(st+1:N,:);
end
B = full(B(1:st,:) - A(1:st,st+1:N)*w1);
vE=ML\B;
timings(1,1)= timings(1,1)+toc(ls);
if nargin < 4 || isempty(maxit), maxit = min(N,50); end;
if nargin < 3 || isempty(tol), tol = 1e-8; end;
if N == 1, fmR = fm(A)*v; err = 0; return; end;
fAb_conv=0;
solver = @(A,b) A\b;    % can use your favorite solver here
zeta = 1*zeros(maxit,1);   % may be optimized
beta = norm(vE);
nu1 = s;
ort=tic;
vEo=vE;
nu1 = s;
[vE,beta]=qr(vE,0);

timings(2,1)=toc(ort);
if adaptive && (strcmp(sh,'conv') || strcmp(sh,'convR') )
    tsh=tic;
    eops.tol=1e-1;
    try
        s1=real(eigs(A,E,1,'SM',eops));
    catch
        s1=-0.1;
    end
    try
        nA=-abs(eigs(A,E,1,'LM',eops));
    catch
        nA=-normest(A);
    end
    timings(4,1)= timings(4,1)+toc(tsh);
end
V = zeros(st,s*2*maxit+1);
AV=V;
V(:,1:nu1) = vE;
ls=tic;
w=MU\vE;
w2=(A(st+1:N,1:st)*w);
switch opts.saveLU_A
    case 1  % use LU factorization
        w1=UA\(LA\w2);
    case 2 % use LU factorization with pivoting
        w1(q,:)=UA\(LA\w2(p,:));
    otherwise %no factor.
        w1=A(st+1:N,st+1:N)\w2(st+1:N,:);
end
w2=A(1:st,1:st)*w - A(1:st,st+1:N)*(w1);
AV(:,1:nu1)=ML\w2;
timings(1,1)= timings(1,1)+toc(ls);


tsh=tic;
R=vE'*AV(:,1:nu1);
theta=eig(R);
timings(4,1)= timings(4,1)+toc(tsh);
H = zeros(2*s*maxit+1,2*s*maxit);
I = speye(N);O=0*Is;
reorth = 1; fmR0 = []; xi = xi(:).';
fmR=zeros(N,1);fmR0a=[];
j=1;jc=1; cplx=0;
curr_start = 1; j = 1; sigma = -1;
while j <= maxit,
    prev_start = curr_start;      % .. first index in the previous step.
    prev_end   = prev_start+nu1-1; % .. last index win the previous step. [curr_start, curr_end]
    real_start = prev_end+1;      % .. first index new columns (for the real part) current step.
    real_end   = real_start+nu1-1; % .. last index new columns (for the real part) current step. [real_start, real_end].
    imag_start = real_end+1;      % .. first index new columns (for the imag. part) current step if shift is complex.
    imag_end   = imag_start+nu1-1; % .. last index new columns(for the imag. part) current step if shift is complex. [imag_start, imag_end].
    
    jms=(j-1)*s+1;j1s=(j+1)*s;js=j*s;js1=js+1;
    if adaptive, %shift computation
        sigma_prev = sigma;
        % look for minimum of nodal rational function sn
        tsh=tic;
        theta(real(theta)>0)=-theta(real(theta)>0);
        if strcmp(sh,'convR'), theta=real(theta); end
        thetao=theta;
        %         theta_aug=[theta;-normest(A)];
        if strcmp(sh,'conv') || strcmp(sh,'convR') %&& length(theta)>2 %
            %           convex hull test set ala
            %           Simoncini/Druskin
            if any(imag(theta)) && length(theta)>2
                theta=[theta;nA];
                ch=convhull(real(theta),imag(theta));
                evH=-theta(ch);
                neH=length(evH); gap=j*s-neH;
                while gap>0, % include enough points from the border
                    neweH=(evH(1:neH-1)+evH(2:neH))/2;
                    evH=[evH;neweH];
                    gap=j*s-length(evH);
                end
                
            else
                evH=sort([-real(theta);-s1;-nA]);
            end
            if j==1, evH=sort(-[s1;nA]); end
            xi_cand=[];
            for k=1:length(evH)-1
                xi_cand=[xi_cand,linspace(evH(k),evH(k+1),500/s)];
            end
        end
        if j==1
            if strcmp(sh,'conv') || strcmp(sh,'convR'),
                gs=-s1*ones(s,1);
            else, gs=inf*ones(s,1);
            end
        else
            gs=kron([xi(2:end)],ones(1,s))';
        end
        sn=ratfunc(xi_cand,thetao,gs);
        [~,jx]=max(abs(sn));
        if real(xi_cand(jx))<0, xi_cand(jx)=-xi_cand(jx); end
        if abs(imag(xi_cand(jx))) / abs( xi_cand(jx))<1e-8
            xi(j+1)=real(xi_cand(jx));
        else
            xi(j+1)=xi_cand(jx);
        end
        timings(4,jc)= timings(4,jc)+toc(tsh);
    end;
    sigma=xi(j+1);
    if( sigma_prev == conj(sigma) )
        % .. Complex conjugate shift, we did this one already, skip
        continue;
    end
    
    ls=tic;
    w=ML*V(:,prev_start:prev_end);
    ww=(A - xi(j+1)*E)\[w;zeros(N-st,size(w,2))];
    w=MU*ww(1:st,:);
    %     w = MU*(solver(A - xi(n+1)*E,ML*V(:,n*s-s+1:n*s)));
    timings(1,jc)= timings(1,jc)+toc(ls);
    ort=tic;
    wR = real(w);
    for it=1:2,
        for kk=1:j
            k1=(kk-1)*nu1+1; k2=kk*nu1;
            gamma=V(:,k1:k2)'*wR;
            H(k1:k2,real_start-nu1:real_end-nu1) = H(k1:k2,real_start-nu1:real_end-nu1)+ gamma;
            wR = wR - V(:,k1:k2)*gamma;
        end
    end
    [V(:,real_start:real_end),H(real_start:real_end, real_start-nu1:real_end-nu1)]=qr(wR,0);

%     [H(real_start:real_end, real_start-nu1:real_end-nu1),V(:,real_start:real_end)]=gram_sh_new(wR);
    %     H(real_start:real_end, real_start-nu1:real_end-nu1)=inv(hinv);
    
    if ~isreal(xi(j+1))
        cplx=1;
        wR=imag(w);
        for it=1:2,
            for kk=1:j+1
                k1=(kk-1)*s+1; k2=kk*s;
                gamma=V(:,k1:k2)'*wR;
                H(k1:k2,imag_start-nu1:imag_end-nu1) = H(k1:k2,imag_start-nu1:imag_end-nu1)+ gamma;
                wR = wR - V(:,k1:k2)*gamma;
            end
        end
        [V(:,imag_start:imag_end),H(imag_start:imag_end, imag_start-nu1:imag_end-nu1),]=qr(wR,0);
        curr_start = real_start;  % .. Where the new columns in Z start
        curr_end   = imag_end;    % .. Where the new columns in Z end
        proj_end   = real_end;    % .. Where the columns used for projection end
    else
        cplx=0;
        % .. The shift is real
        curr_start = real_start;
        curr_end   = real_end;
        proj_end   = prev_end;
    end
    timings(2,jc)= timings(2,jc)+toc(ort);
    ssol=tic;
    dw=curr_end-curr_start+1;
    w(:,1:dw)=MU\V(:,curr_start:curr_end);
    w1(:,1:dw)=A(st+1:N,1:st)*w(:,1:dw);
    switch opts.saveLU_A
        case 1  % use LU factorization
            w1(:,1:dw)=UA\(LA\w1(:,1:dw));
        case 2 % use LU factorization with pivoting
            w1(q,1:dw)=UA\(LA\w1(p,1:dw));
        otherwise %no factor.
            w1(:,1:dw)=A(st+1:N,st+1:N)\w1(:,1:dw);
    end
    w2(:,1:dw)=A(1:st,1:st)*w(:,1:dw) - A(1:st,st+1:N)*w1(:,1:dw);
    AV(:,curr_start:curr_end)=ML\w2(:,1:dw);
    g  = V(:, 1:prev_end)' * AV(:, curr_start:curr_end);
    g3 = V(:, curr_start:curr_end)' * AV(:, 1:curr_end);
    R=[R,g; g3];
    timings(3,jc)=timings(3,jc)+toc(ssol);
    justsolved=0;
    cplxssolve=j>1 && ((xi(j-1)==conj(xi(j))) && mod(j-1,sfreq)==0);
    if (mod(j,sfreq)==0) || cplxssolve
        justsolved=1;
        err(1,jc) = j1s-s;
        if ~fAb_conv
            %%compute small solution
            fmR=expm(R(1:proj_end, 1:proj_end)*max(T))*[beta;zeros(proj_end-nu1,s)];
            df=size(fmR,1)-size(fmR0,1);
            if ~isempty(exact), % rel. change
                err(2,jc) = norm(V(:,1:proj_end)*fmR - exact)/Nexact;
            else %rel. change
                err(2,jc) = norm([fmR0;zeros(df,s)]-fmR)/norm(fmR);
            end;
            if length(T)>1% lower time bound>0
                fmRa=expm(R(1:proj_end, 1:proj_end)*min(T))*[beta;zeros(proj_end-nu1,s)];
                dfa=size(fmRa,1)-size(fmR0a,1);
                if ~isempty(exact), % rel. change
                    err(2,jc) = norm(V(:,1:proj_end)*fmRa - exact)/Nexact;
                else %rel. change
                    err(2,jc) = max(norm([fmR0a;zeros(df,s)]-fmRa)/norm(fmRa),err(2,jc));
                end;
            end
            
            
            fprintf('Step %d -- %e       %d        \n',j,err(2,jc),j*s);
            if err(2,jc) < tol,
                funAb_ss = V(:,1:proj_end)*fmR;
                funAb(pM,:)=ML*funAb_ss;
                if length(T)>1, funAb_ss_a = V(:,1:proj_end)*fmRa;  end
                fAb_conv=1;
            end;
            jc=jc+1;
            fmR0=fmR;
            if length(T)>1,  fmR0a=fmRa; end
        end
        if fAb_conv && tolP>0
            if ~adaptive || (~isnumeric(sh) && strcmp(sh,'im')),
                sh='conv'; if symA && symM, sh='convR'; end
                adaptive=1;
                xi=[1,xi(1:j)];
            end
            if js==size(fmR,1)
                yt=fmR; 
                if length(T)>1, yta=fmRa; end
            else
                yt=V(:,1:proj_end)'*funAb_ss;
                if length(T)>1, yta=V(:,1:proj_end)'*funAb_ss_a; end
            end
            if modG % modified TL Gramian
                eigsops.issym=1; %eigs settings
                eigsops.tol=1e-10;
                eigsops.isreal=1;
                if ~modG_rhs % build the rhs factor of the modified fl CALEs
                    [Db,Eb]=eigs(@(x) [B,funAb]*(flip*([B,funAb]'*x)),st,2*s,'LM',eigsops);
                    e=diag(Eb);
                    [e,idx]=sort(abs(e),'descend');
                    Db=Db(:,idx);
                    beta=(Db(:,1:2*s)*diag(e(1:2*s).^(.5)));
                    switch opts.facM
                        case 1
                            Ba=ML\beta;
                        case 2
                            Ba=ML\beta(pM,:);
                        otherwise
                            Ba=ML\beta;
                    end
                    modG_rhs=1; funAb=beta;
                    if adaptive, sh='conv'; end
                end
                Ba_p=V(:,1:proj_end)'*Ba;
                rhs2=Ba_p*Ba_p';nBf=norm(rhs2);
            else
                if length(T)>1
                    rhs2=yta*yta'-yt*yt';
                else
                    roff=[beta;zeros(proj_end-nu1,s)];
                    rhs2=roff*roff'-yt*yt';
                end
                nBf=norm(rhs2);
            end
            Y = lyap(R(1:proj_end, 1:proj_end),rhs2);
            timings(3,jc)=timings(3,jc)+toc(ssol);
            % computed residual   (exact, in exact arithmetic)
            tres=tic;
            g1 = V(:, 1:curr_end-nu1)' * AV(:, curr_end-nu1+1:curr_end);
            u1 = AV(:, curr_end-nu1+1:curr_end) - V(:, 1:curr_end-nu1) * g1;
            
            g2 = V(:, curr_end-nu1+1:curr_end)' * u1; % p x p matrix
            u2 = u1 - V(:, curr_end-nu1+1:curr_end) * g2;
            
            g3 = qr( u2, 0 ); g3 = triu( g3(1:size(g3, 2), :) ); % p x p matrix; note Y is a proj_end x proj_end matrix, proj_end = k*nu1
            
            HT = H(curr_end-nu1+1:curr_end, proj_end-nu1+1:proj_end); % p x p matrix
            if( proj_end >= 2*nu1 )
                temp = [zeros(nu1, proj_end-nu1) HT] / H(1:proj_end, 1:proj_end); % nu1 x proj_end matrix
                R1 = ( [zeros(nu1, proj_end-2*nu1) -imag(xi(j+1))*HT real(xi(j+1))*HT] / H(1:proj_end, 1:proj_end) - g2*temp ) * Y;
                R2 = -g3*temp*Y;
            else
                temp = [zeros(nu1, proj_end-nu1) HT] / H(1:proj_end, 1:proj_end);
                R1 = ( (real(xi(j+1))*HT) / H(1:proj_end, 1:proj_end) - g2*temp ) * Y;
                R2 = -g3*temp*Y;
            end
            nrmres=norm( [R1; R2] ) / nBf;
            nrmrestot=[nrmrestot,nrmres];
            err(1,jc) = proj_end;
            err(2,jc) = nrmres;
            timings(5,jc)=timings(5,jc)+toc(tres);
            fprintf('Lyap: \t %d\t %e\t %d\n',j,nrmres,j*s)
            
            if (nrmres<tol),
                % factored solution
                if tolsol
                    [uY,sY,~]=svd(Y); [sY,id]=sort(diag(sY),'descend');
                    uY=uY(:,id);
                    is=find(sY>tolsol*sY(1));
                    Y0 = uY(:,is)*diag(sqrt(sY(is)));
                    ls=tic;
                    Z = MU\(V(:,1:size(Y0,1))*Y0);
                    timings(1,jc)= timings(1,jc)+toc(ls);
                else
                    Z=Y;
                end
                V=V(:,1:proj_end);
                break,
                jc=jc+1;
            end
        end
    end;
    if( ~isreal( sigma ) )
        % Two step were done at once.
        j = j + 1;
        curr_start      = imag_start;
        xi(j+1)=conj(xi(j));
    end
    if j >= maxit
        V=V(:,1:proj_end);
        Z=[];
        disp(['Max. iteration number of ' num2str(maxit) ' reached.']);
    end;
    tsh=tic;
    % compute Ritz values of order n
    theta=eig(R(1:proj_end,1:proj_end));
    timings(4,jc)=timings(4,jc)+toc(tsh);
    j=j+1;
end
warning off

function f=ratfunc(x,evH,s)
% evaluate rational function wrt  ritzvalues/shifts
    f=zeros(1,length(x));
    for j=1:length(x)
        f(j)=abs(prod( (x(j)-s)./(x(j)-evH) ));
    end
