function [Z,D,Y,res,niter,timings]=lr_gfadi_r_new(A,B,M,C,F,G,pA,pB,maxiter,rtol)
% function [Z,D,Y,res,niter]=lr_fadi(A,B,F,G,pA,pB,maxiter,rtol,tcrit)
%
% Generate a real low rank matrices Z,Y and D such that X=Z*D*Y' solves the
% Sylvester equation: 
%        
% A*X*C -M*X*B = F*G'                              (1)
%
% The function implements a rearranged version of the low rank fADI method
% as proposed by Benner/Li/Truhar[1] (Algorithm 1).
%
% Inputs: 
% 
% A,B,C,M,F,G      The matrices in the above euqation.
%  pA,pB       a vectors of SUITABLY ORDERED proper shift parameters
% maxiter  maximum iteration number.
% rtol     tolerance for the residual norm based stopping criterion

% Outputs:
%  
% Z,D,Y     The solution factors such that Z*D*Y'=X solves (1).
% res   the vector of residuals in 1st row, steps of measurement in 2nd row
% niter the number of iteration steps taken

%
% [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 Dezember 2012-24

% 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 checked!

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);
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);
nrmZ=0;nrmY=0;
Vt=F; Wt=G;
Z=zeros(n,r*(maxiter+1));Y=zeros(m,r*(maxiter+1));
D=zeros(r*(maxiter+1),r*(maxiter+1));
res=zeros(2,maxiter);
rc=zeros(3,maxiter);rc(1:2,1)=1;
resLan=0;
i=1;j=1;nZo=1;nYo=1;
while i<maxiter+1
  ipA=mod(i+lA-1,lA)+1;
  ipB=mod(i+lB-1,lB)+1;
  ac=pA(ipA);
  bc=pB(ipB);
  ac1=pA(mod(i+lA,lA)+1);
  bc1=pB(mod(i+lB,lB)+1);
  V1=(A-bc*M)\Vt;
  W1=(B-ac*C)'\Wt;
  
  if isreal(ac) && isreal(bc) %case 1: two real shifts
      Z(:,((i-1)*r+1):i*r)=V1;
      Y(:,((i-1)*r+1):i*r)=W1;
      D(((i-1)*r+1):i*r,((i-1)*r+1):i*r)=(bc-ac)*Ir;
      Vt=Vt+(bc-ac)*M*V1; Wt=Wt+(ac-bc)'*C'*W1;
  elseif ~isreal(ac) && ~isreal(bc) %case 2: two complex pairs
      g=real(bc)-real(ac);
      delta=imag(bc)-imag(ac);
      absd=abs(bc-ac)^2;
      Z(:,((i-1)*r+1):(i+1)*r)=[real(V1),imag(V1)];
      Y(:,((i-1)*r+1):(i+1)*r)=[real(W1),imag(W1)];
      D(((i-1)*r+1):(i+1)*r,((i-1)*r+1):(i+1)*r)=[2*g*Ir, (absd+2*imag(ac)*delta)/imag(ac)*Ir;
          (absd-2*imag(bc)*delta)/imag(bc)*Ir,... 
          (g*(2*imag(ac)*imag(bc)+absd))/(imag(ac)*imag(bc))*Ir];
      Vt=Vt+2*g*M*real(V1)+(absd/imag(bc)-2*delta)*M*imag(V1); 
      Wt=Wt-2*g*C'*real(W1)-(absd/imag(ac)+2*delta)*C'*imag(W1);
      bc=conj(bc); ac=conj(ac);
      i=i+1;
  elseif (isreal(ac) && isreal(ac1)) && ~isreal(bc) %case 3a: two real and one complex pair
      g1=real(bc)-ac; g2=real(bc)-ac1;
      delta=ac1^2-2*real(bc)*ac1+abs(bc)^2;
      Y(:,((i-1)*r+1):i*r)=W1;
      Wt = Wt -(g1+g2)*C'*W1;
      Vt = Vt + (g1+g2)*M*real(V1)+...
          (g1*g2-imag(bc)^2)/imag(bc)*M*imag(V1);      
      W1=(B-ac1*C)'\(C'*W1);
      Wt = Wt+delta*C'*W1;
      Z(:,((i-1)*r+1):(i+1)*r)=[real(V1),imag(V1)];
      Y(:,(i*r+1):(i+1)*r)=W1;
      D(((i-1)*r+1):(i+1)*r,((i-1)*r+1):(i+1)*r)=[(g1+g2)*Ir, -delta*Ir;
          (g1*g2-imag(bc)^2)/imag(bc)*Ir,... 
          -g1*delta/imag(bc)*Ir];
      ac=ac1; bc=conj(bc);
      i=i+1;
  elseif (isreal(bc) && isreal(bc1)) && ~isreal(ac) %case 3b: one complex pair and two real 
      g1=bc-real(ac); g2=bc1-real(ac);
      delta=bc1^2-2*real(ac)*bc1+abs(ac)^2; 
      Vt = Vt +(g1+g2)*M*V1;
      Z(:,((i-1)*r+1):i*r)=V1;
      V1 = (A-bc1*M)\(M*V1);
      Vt = Vt +delta*M*V1;
      Wt = Wt -(g1+g2)*C'*real(W1)+...
          (-g1*g2+imag(ac)^2)/imag(ac)*C'*imag(W1);
      Z(:,(i*r+1):(i+1)*r)=V1;
      Y(:,((i-1)*r+1):(i+1)*r)=[real(W1),imag(W1)];
      D(((i-1)*r+1):(i+1)*r,((i-1)*r+1):(i+1)*r)=[(g1+g2)*Ir, (g1*g2-imag(ac)^2)/imag(ac)*Ir;
          delta*Ir,...
          g1*delta/imag(ac)*Ir];
      ac=conj(ac); bc=bc1;
      i=i+1;
  end
  % stopping criteria
  restime=tic;
  [~,Rf]=qr(Vt,0);
  [~,Rg]=qr(Wt,0);
  res(1,j)=norm(Rf*Rg')/res0;
  timings(2,i)=toc(restime);
  res(2,j) = i; j=j+1;
 
  fprintf(1,['step: %4d  n.res.: %e \n'],i,res(1,j-1));     
  if res(1,j-1)<rtol
    fprintf('\n\n');
    break;
  end
  i=i+1;
end
niter=i;
res=res(:,1:j-1);