function [Z,D,Y,res,niter,timings]=lr_adi_sylv_ex_cpx(A,B,M,C,F,G,pA,pB,maxiter,rtol)
% function [Z,D,Y,res,niter]=lr_adi_sylv_ex_cpx(A,B,M,C,F,G,pA,pB,maxiter,rtol,tcrit)
%
% Generate a low rank matrices Z,Y and diagonal D such that X=Z*D*Y' solves the
% Sylvester equation:
%
% A*X*C + E*X*B + F*G'=0                              (1)
%
% The function implements the generalized low rank Sylvester-ADI method,
% see [1,2,3] (also called (g-)fADI).
%
% Inputs:
%
% A,B,E,C,F,G      The matrices in the above equation.
%  pA,pB       vectors of shift parameters
% maxiter  maximum iteration number.
% rtol     tolerance for the residual norm based stopping criterion

% Outputs:
%
% Z,D,Y     The solution factors Z,D,Y such that Z*D*Y'=X solves (1).
% res       the vector of residuals
% niter     the number of iteration steps taken
% timings   comp.times of various stages
%               (1st row: full step, 2nd row: res.norm computation)
%
% [1] Benner/Li/Truhar, On the ADI method for Sylvester equations, J. Comput. Appl. Math.,
%                       346 233(4):1035–1045, 2009.
% [2] Benner/Kuerschner, Computing Real Low-rank Solutions of Sylvester equations by the Factored ADI
%                       Method, Comput. Math. Appl., 67(9):1656–1672, 2014.
% [3] Kuerschner, Inexact linear solves in the low-rank ADI iteration for
%      large Sylvester equations, Arxiv e-print 2312.02891,  2023

% Patrick Kuerschner, 2016-2024

% 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 parameters not fully checked!

% various sizes & dimensions
n=size(A,1);
m=size(B,1);
r=size(F,2);
lA=length(pA);
lB=length(pB);
In=speye(n);
Im=speye(m);
Ir=speye(r);
% U=zeros;
tadi=tic;
%starting residual norm
[~,Rg]=qr(F,0);
[~,Rf]=qr(G,0);
res0=norm(Rf*Rg');
timings(1,1)=toc(tadi);

ac=pA(1);
bc=pB(1);
Vt=F; Wt=G;
Z=[];Y=[];D=[];
res=zeros(2,maxiter);

opts.isreal=false;
opts.issym=true;
timings=zeros(2,maxiter);
timings(1,1)=toc(tadi);
for i=1:maxiter %Sylvester-ADI Loop
    tadi=tic;
    %select pair of shifts
    ipA=mod(i+lA-1,lA)+1;
    ipB=mod(i+lB-1,lB)+1;
    ap=ac;
    ac=pA(ipA);
    bp=bc;
    bc=pB(ipB);
    % solve lin.sys.
    V1=(A+bc*M)\Vt;
    W1=(B+ac*C)'\Wt;
    %augment low-rank factor
    Z=[Z V1];
    Y=[Y W1];
    gam=-(bc+ac);
    D=blkdiag(D,gam*Ir);
    %update residual factors
    gEV=gam*MV(M,V1); gCW=conj(gam)*MV(C',W1);
    % residual factors
    Vt=Vt+gEV; Wt=Wt+gCW;
    % residual norm


    restime=tic;
    [~,Rf]=qr(Vt,0);
    [~,Rg]=qr(Wt,0);
    res(1,i)=norm(Rf*Rg')/res0;
    timings(2,i)=toc(restime);

    fprintf(1,['step: %4d  n.resLR: %2.4e \n'],i,res(1,i));

    %are we done?
    if res(1,i)<rtol
        fprintf('\n\n');
        timings(1,i)=timings(1,i)+toc(tadi);
        break;
    end
    timings(1,i)=timings(1,i)+toc(tadi);
end
niter=i;
timings=timings(:,1:niter);
res=res(:,1:niter);