function [kkt1,kkt2,lambda_1opt,fval,xopt,time,LNGM,xopt1]=Algorithm2(A,a,b,beta,delta,n,lambdamin2)
% This function soves the following eTRS by Algorithm 2 (The eTRS does not eonjoy strong duality)

%  min x'Ax+2a'x
%      x'x<=delta
%      b'x<=beta

% Outputs:
% kkt1 and kkt2: first order stationary conditions
% lambda_1opt: optimal lagrange multiplier
% fval: optimal value of eTRS
% xopt: optimal solution of eTRS
% time: time taken for Algorithm 2
% LNGM:  LNGM=1 if optimal solution of eTRS is the LNGM of corresponding TRS otherwise equals 0

global flag_con  bbar bl W change xhat lambdamin
%============================================================
%==  computing local-nonglobal minimizer of TRS (LNGM)
LNGM=0;
tic
fla=0;  % fla will changed to one if there is no LNGM for TRS
flag_con=0;
M1=[sparse(n,n) speye(n);speye(n) sparse(n,n)];
opts.v0=randn(2*n,1);
opts.issym=1;
opts.isreal=1;
[v,d]=eigs(@M0,2*n,-M1,2,'lr',opts);
if (abs(imag(d(2,2))))>1e-8
    fprintf('Second eigenvalue of M(lambda) is complex \n')
    if real(d(2,2))<=-lambdamin2
        fla=1;
        time1=toc;
        fprintf('fla=1: there is no local nonglobal minimizer\n')
    else
        [~,d,flagsecon]=realM0n(v(:,1),d(1,1),lambdamin2);
        if flagsecon==0,
            fla=1;
            time1=toc;
            fprintf('fla=1: there is no local nonglonal minimizer \n');
        elseif (d>=0)&&((d+lambdamin2)>=1e-8)&&((d+lambdamin1)<-1e-8)
            df=d;
            opts.v0=randn(2*n,1);
            opts.issym=0;
            opts.isreal=1;
            [v,d]=eigs(@Asigma,2*n,1,df,opts);
            xopt1=-sign(a'*v(n+1:end))*sqrt(delta)*v(1:n)/norm(v(1:n));
            fval1=xopt1'*A*xopt1+2*a'*xopt1;
            time1=toc;
            fprintf('eigenvector for the second real eige of M(lambda) is computed\n')
            lambda_1LNGM=d;
            kkt1=norm((A+lambda_1LNGM*speye(n))*xopt1+a,inf);
            kkt2=lambda_1LNGM*(norm(xopt1)^2-delta);
            fprintf('kkt1 and kkt2 for LNGM are: %g  %g \n',kkt1,kkt2);
        else
            fla=1;
            time1=toc;
            fprintf('fla=1: there is no local nonglonal minimizer \n');
            
        end
    end
elseif (d(2,2)>1e-12)&&((d(2,2)+lambdamin2)>=1e-8)
    xopt1=-sign(a'*v(n+1:end,2))*sqrt(delta)*v(1:n,2)/norm(v(1:n,2));
    fval1=xopt1'*A*xopt1+2*a'*xopt1;
    time1=toc;
    lambda_1LNGM=d(end,end);
    kkt1=norm((A+lambda_1LNGM*speye(n))*xopt1+a,inf);
    kkt2=lambda_1LNGM*(norm(xopt1)^2-delta);
    fprintf('kkt1 and kkt2 for LNGM are: %g  %g \n',kkt1,kkt2);
elseif (abs(d(2,2))<=1e-12)&&((d(2,2)+lambdamin2)>=1e-8)
    xopt1=-sign(a'*v(n+1:end,2))*sqrt(delta)*v(1:n,2)/norm(v(1:n,2));
    fval1=xopt1'*A*xopt1+2*a'*xopt1;
    time1=toc;
    lambda_1LNGM=d(end,end);
    kkt1=norm((A+lambda_1LNGM*speye(n))*xopt1+a,inf);
    kkt2=lambda_1LNGM*(norm(xopt1)^2-delta);
    fprintf('kkt1 and kkt2 for LNGM are: %g  %g \n',kkt1,kkt2);
else
    time1=toc;
    fprintf('fla=1: there is no  local minimzer\n')
    fla=1;
end
%======================================================================
if fla==1
    xopt1=zeros(n,1);
end

if ((xopt1'*b-beta)>1e-6)||(fla==1)
    % LNGM is not feasible for lin con or there is no local minima,  then optimal solution of eTRS is
    % the solution of projected TRS
    %=============================================================
    %=== solving the projected TRS problem =======================
    tic
    flag_con=1;          % indicates we are solving the projected eTRS problem (eTRS with active linear constraint)
    [bl,ind]=max(abs(b));
    ind=ind(1);
    if bl<1e-4
        fprintf('Vector b is almost zero\n');
    else
        bl=b(ind);        % bl : |bl|>=|b_i|  , i=1,...,n
        bbar=b;
        bbar(ind)=[];        % bbar is obtaind by removing bl from b
        W=[bl*speye(ind-1) sparse(ind-1,n-ind);-bbar'; sparse(n-ind,ind-1) bl*speye(n-ind)]; %  W is a basis for null(b)
        if beta==0
            % xhat=0   is a feasible for b'x=beta
            opts.v0=randn(2*(n-1),1);
            opts.issym=1;
            opts.isreal=1;
            M1=[sparse(n-1,n-1) speye(n-1);speye(n-1) sparse(n-1,n-1)];
            [v,d]=eigs(@M0,2*(n-1),-M1,1,'lr',opts);
            bnorm=norm(bbar)^2;
            coef=(1+sqrt((bl^2)/(bl^2+bnorm)))/bnorm;
            v1=(1/abs(bl))*(v(1:n-1)-coef*bbar*(bbar'*v(1:n-1)));
            v2=(1/abs(bl))*(v(n:end)-coef*bbar*(bbar'*v(n:end)));
            v=[v1;v2];
            if norm(v(1:n-1))<=1e-6   %1e-4 does not work
                fprintf('The projected TRS problem is a hard case instance\n')
                lambdamin=real(d);
                v=v(n:end)/norm(v(n:end));
                [base,alpha]=basisnull2(v,lambdamin);
                base=B_orthog(base);
                q=pcg(@Hfunc,-W'*a,1e-14,1000);
                Bv=W'*(W*v);
                vBv=v'*Bv;
                xyprod = Bv'*q;  % normalized inner product
                ndiff = delta;
                % eta=roots([vBv 2*xyprod q'*(W'*(W*q))-xhat'*(W*change)-delta+norm(xhat)^2]);
                % tau=eta(1);
                sxy = sign(xyprod);
                if sxy == 0,
                    sxy = 1;
                end
                tau = ndiff/(xyprod + sxy*sqrt(xyprod^2 + vBv*ndiff));
                xopt0 = q + tau*v;
                xopt=W*xopt0;
                time2=toc;
                opts.v0=randn(n-1,1);
                opts.issym=1;
                opts.isreal=1;
                [~,eigc]=eigs(@chek,n-1,1,'sa',opts);
                fprintf(' smallest eigenvalue of (A+lambda_g B)=%g\n',eigc);
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))+W'*a,inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g  %g \n',kkt1,kkt2);
            else
                normBy1=sqrt(v(1:n-1)'*(W'*(W*v(n:end))));
                xopt0=-sign(a'*(W*v(n:end)))*sqrt(delta)*v(1:n-1)/normBy1;
                xopt=W*xopt0;
                time2=toc;
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))+W'*a,inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g  %g \n',kkt1,kkt2);
            end
            
        else
            bnorm=norm(bbar)^2;
            xhat=(beta/(norm(b)^2))*b;    %xhat based on formula in page 12
            % change of variable to transform nTRS to homogenous one
            change=(1/(bl^2))*(W'*xhat-(1/(bl^2+bnorm))*bbar*(bbar'*(W'*xhat)));
            opts.v0=randn(2*(n-1),1);
            opts.issym=1;
            opts.isreal=1;
            M1=[sparse(n-1,n-1) speye(n-1);speye(n-1) sparse(n-1,n-1)];
            [v,d]=eigs(@M0,2*(n-1),-M1,1,'lr',opts);
            coef=(1+sqrt((bl^2)/(bl^2+bnorm)))/bnorm;
            v1=(1/abs(bl))*(v(1:n-1)-coef*bbar*(bbar'*v(1:n-1)));
            v2=(1/abs(bl))*(v(n:end)-coef*bbar*(bbar'*v(n:end)));
            v=[v1;v2];
            if norm(v(1:n-1))<=1e-6 %1e-4 does not work
                fprintf('The projectd TRS problem is a hard case instance\n')
                lambdamin=real(d);
                v=v(n:end)/norm(v(n:end));
                [base,alpha]=basisnull2(v,lambdamin);
                base=B_orthog(base);
                g=-(W'*(A*(W*change)))+W'*((A*xhat)+a);
                q=pcg(@Hfunc,-g,1e-14,1000);
                Bv=W'*(W*v);
                vBv=v'*Bv;
                xyprod = Bv'*q;  % normalized inner product
                ndiff = xhat'*(W*change)+delta-norm(xhat)^2;
                %             eta=roots([vBv 2*xyprod q'*(W'*(W*q))-xhat'*(W*change)-delta+norm(xhat)^2]);
                %             tau=eta(1);
                sxy = sign(xyprod);
                if sxy == 0,
                    sxy = 1;
                end
                tau = ndiff/(xyprod + sxy*sqrt(xyprod^2 + vBv*ndiff));
                xopt0 = q + tau*v;
                time2=toc;
                opts.v0=randn(n-1,1);
                opts.issym=1;
                opts.isreal=1;
                [~,eigc]=eigs(@chek,n-1,1,'sa',opts);
                fprintf(' smallest eigenvalue of (A+lambda_g B)=%g\n',eigc);
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))-(W'*(A*(W*(change))))+(W'*(A*xhat+a)),inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-xhat'*(W*change)+norm(xhat)^2-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g %g\n',kkt1,kkt2);
                
                
            else
                gy2=(xhat'*A+a')*(W*(v(n:end)))-change'*(W'*(A*(W*v(n:end))));
                normBy1=sqrt(v(1:n-1)'*(W'*(W*(v(1:n-1)))));
                gamma=sqrt(xhat'*(W*change)+delta-norm(xhat)^2);
                xopt0=-sign(gy2)*gamma*v(1:n-1)/normBy1;
                time2=toc;
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))-(W'*(A*(W*(change))))+(W'*(A*xhat+a)),inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-xhat'*(W*change)+norm(xhat)^2-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g %g \n',kkt1,kkt2);
            end
            xopt0=xopt0-change;
            xopt=xhat+W*xopt0;
        end
    end
    fval=xopt'*A*xopt+2*a'*xopt;
    fval=real(fval);
    time=time1+time2;
else
    %=============================================================
    %=== solving the projected TRS problem =======================
    tic
    flag_con=1;          % indicates we are solving the projected eTRS problem (eTRS with active linear constraint)
    [bl,ind]=max(abs(b));
    ind=ind(1);
    if bl<1e-4
        fprintf('Vector b is almost zero\n');
    else
        bl=b(ind);        % bl : |bl|>=|b_i|  , i=1,...,n
        bbar=b;
        bbar(ind)=[];        % bbar is obtaind by removing bl from b
        W=[bl*speye(ind-1) sparse(ind-1,n-ind);-bbar'; sparse(n-ind,ind-1) bl*speye(n-ind)]; %  W is a basis for null(b)
        if beta==0
            % xhat=0   is a feasible for b'x=beta
            opts.v0=randn(2*(n-1),1);
            opts.issym=1;
            opts.isreal=1;
            M1=[sparse(n-1,n-1) speye(n-1);speye(n-1) sparse(n-1,n-1)];
            [v,d]=eigs(@M0,2*(n-1),-M1,1,'lr',opts);
            bnorm=norm(bbar)^2;
            coef=(1+sqrt((bl^2)/(bl^2+bnorm)))/bnorm;
            v1=(1/abs(bl))*(v(1:n-1)-coef*bbar*(bbar'*v(1:n-1)));
            v2=(1/abs(bl))*(v(n:end)-coef*bbar*(bbar'*v(n:end)));
            v=[v1;v2];
            if norm(v(1:n-1))<=1e-6  %1e-4 does not work
                fprintf('The projected TRS problem is a hard case instance\n')
                lambdamin=real(d);
                v=v(n:end)/norm(v(n:end));
                [base,alpha]=basisnull2(v,lambdamin);
                base=B_orthog(base);
                q=pcg(@Hfunc,-W'*a,1e-14,1000);
                Bv=W'*(W*v);
                vBv=v'*Bv;
                xyprod = Bv'*q;  % normalized inner product
                ndiff = delta;
                % eta=roots([vBv 2*xyprod q'*(W'*(W*q))-xhat'*(W*change)-delta+norm(xhat)^2]);
                % tau=eta(1);
                sxy = sign(xyprod);
                if sxy == 0,
                    sxy = 1;
                end
                tau = ndiff/(xyprod + sxy*sqrt(xyprod^2 + vBv*ndiff));
                xopt0 = q + tau*v;
                xopt=W*xopt0;
                time2=toc;
                opts.v0=randn(n-1,1);
                opts.issym=1;
                opts.isreal=1;
                [~,eigc]=eigs(@chek,n-1,1,'sa',opts);
                fprintf(' smallest eigenvalue of (A+lambda_g B)=%g\n',eigc);
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))+W'*a,inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g  %g \n',kkt1,kkt2);
            else
                normBy1=sqrt(v(1:n-1)'*(W'*(W*v(n:end))));
                xopt0=-sign(a'*(W*v(n:end)))*sqrt(delta)*v(1:n-1)/normBy1;
                xopt=W*xopt0;
                time2=toc;
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))+W'*a,inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g  %g \n',kkt1,kkt2);
            end
            
        else
            bnorm=norm(bbar)^2;
            xhat=(beta/(norm(b)^2))*b;    %xhat based on formula in page 12
            % change of variable to transform nTRS to homogenous one
            change=(1/(bl^2))*(W'*xhat-(1/(bl^2+bnorm))*bbar*(bbar'*(W'*xhat)));
            opts.v0=randn(2*(n-1),1);
            opts.issym=1;
            opts.isreal=1;
            M1=[sparse(n-1,n-1) speye(n-1);speye(n-1) sparse(n-1,n-1)];
            [v,d]=eigs(@M0,2*(n-1),-M1,1,'lr',opts);
            coef=(1+sqrt((bl^2)/(bl^2+bnorm)))/bnorm;
            v1=(1/abs(bl))*(v(1:n-1)-coef*bbar*(bbar'*v(1:n-1)));
            v2=(1/abs(bl))*(v(n:end)-coef*bbar*(bbar'*v(n:end)));
            v=[v1;v2];
            if norm(v(1:n-1))<=1e-6  %1e-4 does not work
                fprintf('The projectd TRS problem is a hard case instance\n')
                lambdamin=real(d);
                v=v(n:end)/norm(v(n:end));
                [base,alpha]=basisnull2(v,lambdamin);
                base=B_orthog(base);
                g=-(W'*(A*(W*change)))+W'*((A*xhat)+a);
                q=pcg(@Hfunc,-g,1e-14,1000);
                Bv=W'*(W*v);
                vBv=v'*Bv;
                xyprod = Bv'*q;  % normalized inner product
                ndiff = xhat'*(W*change)+delta-norm(xhat)^2;
                %             eta=roots([vBv 2*xyprod q'*(W'*(W*q))-xhat'*(W*change)-delta+norm(xhat)^2]);
                %             tau=eta(1);
                sxy = sign(xyprod);
                if sxy == 0,
                    sxy = 1;
                end
                tau = ndiff/(xyprod + sxy*sqrt(xyprod^2 + vBv*ndiff));
                xopt0 = q + tau*v;
                time2=toc;
                opts.v0=randn(n-1,1);
                opts.issym=1;
                opts.isreal=1;
                [~,eigc]=eigs(@chek,n-1,1,'sa',opts);
                fprintf(' smallest eigenvalue of (A+lambda_g B)=%g\n',eigc);
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))-(W'*(A*(W*(change))))+(W'*(A*xhat+a)),inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-xhat'*(W*change)+norm(xhat)^2-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g %g\n',kkt1,kkt2);
                
                
            else
                gy2=(xhat'*A+a')*(W*(v(n:end)))-change'*(W'*(A*(W*v(n:end))));
                normBy1=sqrt(v(1:n-1)'*(W'*(W*(v(1:n-1)))));
                gamma=sqrt(xhat'*(W*change)+delta-norm(xhat)^2);
                xopt0=-sign(gy2)*gamma*v(1:n-1)/normBy1;
                time2=toc;
                lambda_1opt=d;
                kkt1=norm(W'*(A*(W*(xopt0)))+d*(W'*(W*(xopt0)))-(W'*(A*(W*(change))))+(W'*(A*xhat+a)),inf);
                kkt2=lambda_1opt*(xopt0'*(W'*(W*(xopt0)))-xhat'*(W*change)+norm(xhat)^2-delta);
                fprintf('kkt1 and kkt2 for solution of projected TRS are: %g %g \n',kkt1,kkt2);
            end
            xopt0=xopt0-change;
            xopt=xhat+W*xopt0;
        end
    end
    fval2=xopt'*A*xopt+2*a'*xopt;
    if fla==0
        tic
        [~,i]=min([real(fval1), real(fval2)]);
        if i==1
            xopt=xopt1;
            LNGM=1;     %indicates the global solution of eTRS is the local minimzer
            time3=toc;
            time=time1+time2+time3;
            fprintf('Main algorithm indicates LNGM is the global found')
            lambda_1opt=lambda_1LNGM;
            kkt1=norm((A+lambda_1opt*speye(n))*xopt+a,inf);
            kkt2=lambda_1opt*(norm(xopt)^2-delta);
        else
            time=time1+time2;
        end
    else
        time=time1+time2;
    end
    fval=xopt'*A*xopt+2*a'*xopt;
    fval=real(fval);
end
    %======================================================================
    function [v] = Hfunc(u)
    
    if flag_con==0
        v=A*u+lambdamin*u+(base*(alpha.*(base'*u)));
    elseif flag_con==1
        v=W'*(A*(W*u))+lambdamin*(W'*(W*u))+(base*(alpha.*(base'*u)));
    end
    end
end
