function [kkt1,kkt2,lambda_1opt,strong,LNGM,time,timestrong,fval,xopt,xopt1]=Algorithm3(A,a,b,beta,delta,n)

% This function solves the following problem (eTRS) by Algorithm 3 (efficiently solving a general eTRS; exploiting
%  strong duality )

%      min x'Ax+2a'x
%          x'x<=delta   (eTRS)
%          b'x<=beta
% Outputs:
% kkt1 and kkt2: first order stationary conditions
% lambda_1opt: optimal lagrange multiplier
% strong: strong=1 if eTRS enjoyes strong duality otherwise equals 0 
% LNGM:  LNGM=1 if optimal solution of eTRS is the LNGM of corresponding TRS otherwise equals 0
% time: time taken for Algorithm 3 including the time taken for verifying strong duality property
% timestrong: time taken for verifying strong duality property
% fval: optimal value of eTRS
% xopt: optimal solution of eTRS
warning off
%===============================================================
LNGM=0;    % LNGM will changed to one if the global solution of eTRS is the local nonglobal minimizer of corresponding TRS
timestrong=0;
%===============================================================
tic
opts.v0=full(sum(A))';
opts.issym=1;
opts.isreal=1;
[vh,lambdamin1]=eigs(A,1,'sa',opts);   % computes minimum eigenvalue of A (lambdamin1)
vh=vh/norm(vh);
if lambdamin1>=0
    t=toc;
    timestrong=timestrong+t;
    strong=1;
    fprintf( 'A is positive (semi) definite ');
    [kkt1,kkt2,lambda_1opt,fval,xopt,time_alg1]=Algorithm1(A,a,b,beta,delta,n);
    time=time_alg1+timestrong;
else
    [base,alpha,lambdamin2]=base_eig(A,lambdamin1,vh);  %lambdamin2 is the second smallest eigenvalue of A
    if size(base,2)>=2
        t=toc;
        timestrong=timestrong+t;
        strong=1;
        [kkt1,kkt2,lambda_1opt,fval,xopt,time_alg1]=Algorithm1(A,a,b,beta,delta,n,base,alpha);
        time=time_alg1+timestrong;
    elseif (abs(vh'*b))<1e-8
        t=toc;
        strong=1;
        timestrong=timestrong+t;
        [kkt1,kkt2,lambda_1opt,fval,xopt,time_alg1]=Algorithm1(A,a,b,beta,delta,n,base,alpha);
        time=time_alg1+timestrong;
    elseif ((-vh'*a)/(vh'*b))<=1e-8
        t=toc;
        strong=1;
        timestrong=timestrong+t;
        [kkt1,kkt2,lambda_1opt,fval,xopt,time_alg1]=Algorithm1(A,a,b,beta,delta,n,base,alpha);
        time=time_alg1+timestrong;
    else
        mu=(-vh'*a)/(vh'*b);
        hand=-a-mu*b;
        [qp,~]=pcg(@private,hand,1e-14,100000);
%         fprintf('pcg from MainAlgorithm time %g \n',toc)
        step=roots([1 2*qp'*vh qp'*qp-delta]);
        x1=qp+step(1)*vh;
        x2=qp+step(2)*vh;
        if (abs((b'*x1-beta)*(b'*x2-beta)))<1e-10
            t=toc;
            timestrong=timestrong+t;
            strong=1;
            if (abs(b'*x1-beta))<1e-8
                xopt=x1;
                time=timestrong;
            else
                xopt=x2;
                time=timestrong;
            end
        elseif ((b'*x1-beta)*(b'*x2-beta))>=1e-10
            t=toc;
            timestrong=timestrong+t;
            strong=1;
            [kkt1,kkt2,lambda_1opt,fval,xopt,time_alg1]=Algorithm1(A,a,b,beta,delta,n,base,alpha);
            time=time_alg1+timestrong;
        else
            t=toc;
            timestrong=timestrong+t;
            strong=0;
            [kkt1,kkt2,lambda_1opt,fval,xopt,time_alg2,LNGM,xopt1]=Algorithm2(A,a,b,beta,delta,n,lambdamin2);
            time=timestrong+time_alg2;
        end
    end
end
%=============================================================
    function f=private(x)
        nA = normest(A,1);
        nAn=nA/sqrt(n);
        f=A*x-lambdamin1*x+nAn*(vh*(vh'*x));
    end
end
