function f=M0(x)
% This function computes M0*x
global bbar bl W a flag_con A delta n beta change xhat
%===============================================================
if flag_con==0
    f=[-x(1:n)+A*x(n+1:end);A*x(1:n)-a*(a'*x(n+1:end))/delta];
    
elseif (flag_con==1)&(beta==0)
    %xhat=0
    Binsx2=Binvsqrt(x(n:end));
    part12=Binvsqrt(W'*(A*(W*Binsx2)));
    Binsx1=Binvsqrt(x(1:n-1));
    part21=Binvsqrt(W'*(A*(W*Binsx1)));
    abar=W'*a;
    Binsabar=Binvsqrt(abar);
    part22=-((abar'*Binsx2)/delta)*Binsabar;
    f=[-x(1:n-1)+part12;part21+part22];
else
    Binsx2=Binvsqrt(x(n:end));
    part12=Binvsqrt(W'*(A*(W*Binsx2)));
    Binsx1=Binvsqrt(x(1:n-1));
    part21=Binvsqrt(W'*(A*(W*Binsx1)));
    abar=-(W'*(A*(W*change)))+(W'*((A*xhat)+a));
    Binsabar=Binvsqrt(abar);
    del=xhat'*(W*change)+delta-norm(xhat)^2;
    part22=-((abar'*Binsx2)/del)*Binsabar;
    f=[-x(1:n-1)+part12;part21+part22];
    
end

    function f=Binvsqrt(x)
        % This function computes B^-{1/2}*x which B=W'*W=bl^2*eye(n-1)+bbar*bbar'
        %===============================================================
        bnorm=norm(bbar)^2;
        coef=(1+sqrt((bl^2)/(bl^2+bnorm)))/bnorm;
        f=(1/abs(bl))*(x-coef*bbar*(bbar'*x));
    end
end

