  function  output=MSFprecondfun(mode,n2,npmmpmc,dxy,Qy,diagprec)
% function  output=MSFprecondfun(mode,n2,npmmpmc,dxy,Qy,diagprec)
%
%  INPUT:
%   dxy = ????? contains the vector to be multiplied in CG method
%   
global indsE indsEc X Z;


n=sqrt(n2);   % number of nodes
mm=length(indsE);
mc=(n*(n-1)/2)-mm;    % complement of number of arcs/edges


if mode==2,
    dM=reshape(dxy,n,n);
     %%%%%%%%zops operator
    Zs=Z*dM;    % adjoint multiplication by Z
    Zs=(Zs+Zs')/2;    % adjoint multiplication by Z
    zopsout=diag(Zs); %% adjoint of partial diagonal operator
    zopsout=zopsout(2:n)-zopsout(1);   % vector length n-1
    zopsout=[zopsout
             sqrt(2)*Zs(indsEc) ];  % final vector length n-1 + mc
     %%%%%%%%xops operator
    Xs=dM*X;
    Xs=(Xs+Xs')/2;
    xopsout=[ trace(Xs)
              sqrt(2)*Xs(indsE)];

    output=[zopsout
            xopsout ];
    output=output./diagprec;
%%%%%%%%%%%%%%%%%%
else          % mode ==1
   dxy=dxy./diagprec;
   dd=dxy(1:n-1);
   dx=dxy(n:n+mc-1);
   dz=dxy(n+mc);
   dy=dxy(n+mc+1:end);
%%%%Zopout is zop, Z operator
   Zopout=zeros(n);
   Zopout(indsEc)=dx/sqrt(2);
   Zopout=Zopout+Zopout';
   Zopout(1:n+1:end)=[-sum(dd);dd];
   Zopout=Z*Zopout;
%%%%Xopout is xop, X operator
   Xopout=zeros(n);
   Xopout(indsE)=dy/sqrt(2);
   Xopout=Xopout+Xopout';
   Xopout(1:n+1:end)=diag(Xopout)+dz;
   Xopout=Xopout*X;
   output=Zopout+Xopout;
   output=output(:);
end
