function r = syl_r_norm(A,B,F,G,Z,D,Y,m,M,C,W,T)
% this function estimates the 2-norm of the  nxm Sylvester residual
% matrix R=AXC+MCB+fg' w.r.t. a low-rank X=ZDY'
% by a call to matlab's eigs() routine using MatVecs with R'R
% 
% Input: A,B,E,C,F,G   Coefficient matrices in the above euqation.
%           Z,D,Y          factors of approx. sol  
%         m             col.-dim. of R
%       W,T             for initial vector calc. (experimental)

% 
% Copyright (C) 2016-2024 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

if isempty(M) && isempty(C) % M, C are identity matrices?
    R = @(x) A*(Z*(D*(Y'*x)))+Z*(D*(Y'*(B*x)))+F*(G'*x);
    Rt = @(x) B'*(Y*(D'*(Z'*x)))+Y*(D'*(Z'*(A'*x)))+G*(F'*x);
else
    R = @(x) A*(Z*(D*(Y'*(C*x))))+M*(Z*(D*(Y'*(B*x))))+F*(G'*x);
    Rt = @(x) +B'*(Y*(D'*(Z'*(M'*x))))+C'*(Y*(D'*(Z'*(A'*x))))+G*(F'*x);
end
if nargin>10 && (~isempty(W) && ~isempty(T))
    [u,s]=eigs(@(x)(W'*W)*(T'*T)*x,size(T,2),1,'LM');
    opts.v0=T*u;
else
    opts.v0=ones(m,1);
end    
opts.isreal=(isreal(Z) && isreal(Y) && isreal(D)); %real data?
opts.issym=true; % R'R is symmetric / hermitian
opts.tol=1e-2;
opts.maxit=5;
r=sqrt(abs(eigs(@(x) Rt(R(x)),m,1,'LM',opts))); % call to eigs

