function [kkt1,kkt2,lambda_1opt,fval,xopt,time]=Algorithm1(A,a,b,beta,delta,n,base,alpha)

% This function solves the following eTRS by Algorithm 1(solving eTRS under strong duality property)
%  min x'Ax+2a'x
%      x'x<=delta
%      b'x<=beta

% Outputs:
% kkt1 and kkt2: first order stationary conditions
% lambda_1opt: optimal lagrange multiplier 
% time: time taken for Algorithm 1
% fval: optimal value of eTRS
% xopt: optimal solution of eTRS

global flag_con  bbar bl W change xhat lambdamin
%======================================================================
tic
if nargin<=6
    ps=1;
else
    ps=0;
end
i=0;
flag_con=0;            % indicates we are solving TRS by Algorithm 4
cand=0;                % cand will changed to one if x0 from step 1 of Algorithm 4 is feasible for TRS
if ps==1
    [x0,~]=pcg(A,-a,1e-14,1000);
    if (norm(x0)^2)<=delta
        cand=1;
        fval0=x0'*A*x0+2*x0'*a;
    end
    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,1,'lr',opts);
    if norm(v(1:n))<=1e-4
        fprintf(' TRS is a hard case instance\n');
        lambdamin=real(d);
        v=v(n+1:end)/norm(v(n+1:end));
        [base,alpha,~]=base_eig(A,-lambdamin,v);
        [q,~]=pcg(@Hfunc,-a,1e-14,1000);
        qv=q'*v;
        tu=(delta-(norm(q))^2)/(qv+(sign(qv))^2*(qv^2+delta-norm(q)^2)^(1/2));
        xopt=q+tu*v;
        fval1=xopt'*A*xopt+2*a'*xopt;
    else
        xopt=-sign(a'*v(n+1:end))*sqrt(delta)*v(1:n)/norm(v(1:n));
        fval1=xopt'*A*xopt+2*a'*xopt;
    end
    if cand==1
        [~,i]=min(real([fval0 fval1]));
        if i==1
            xopt=x0;
        end
    end
    
else
    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,1,'lr',opts);
    if norm(v(1:n))<=1e-4
        fprintf(' TRS is a hard case instance\n')
        lambdamin=real(d);
        [q,~]=pcg(@Hfunc,-a,1e-14,1000);
        qv=q'*v(n+1:end);
        tu=(delta-(norm(q))^2)/(qv+(sign(qv))^2*(qv^2+delta-norm(q)^2)^(1/2));
        xopt=q+tu*v(n+1:end);
    else
        xopt=-sign(a'*v(n+1:end))*sqrt(delta)*v(1:n)/norm(v(1:n));
    end
end
time1=toc;
if i==1
    lambda_1opt=0;
else
    lambda_1opt=d;
end
kkt1=norm((A+lambda_1opt*speye(n))*xopt+a,inf);
kkt2=lambda_1opt*(norm(xopt)^2-delta);
fval=xopt'*A*xopt+2*a'*xopt;
%==================================================================
% if optimal solution of TRS (xopt) is feasible for b'*x<=beta, then it is also optimal for eTRS.
time2=0;
tic
if (xopt'*b-beta)>1e-6
    %---------------------------------------------------------------
    %--------------- solving the projected TRS problem -----------------
    flag_con=1;          % indicates we are solving the projected TRS 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 eTRS 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 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);
            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 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 eTRS 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 eTRS are: %g %g \n',kkt1,kkt2);
            end
            xopt0=xopt0-change;
            xopt=xhat+W*xopt0;
        end
    end
    fval=xopt'*A*xopt+2*a'*xopt;
end
time=time1+time2;

%======================================================================
    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

