%Federica Tinti (ftinti@dm.unipd.it)
%Department of Mathematics of Ferrara (ITALY)
%
%Extragradient method (Marcotte's versions)
%
clear all
stringa=input('name of input test problem=', 's');
eval(stringa);
%parameters
P=eye(n);
eps=1e-4;
rx=0;
alpha=1;
beta=0.7;
%
iter=0;      
nf=1;               
np=0;				
%
tempo=cputime;
fx=feval(fs,x);
while norm(rx)>eps | iter==0 | flag==1
        %first projection
        f=-(x-alpha*fx);
        options=optimset('LargeScale','off','display','off');
        [xx]=quadprog(P,f,A,b,Aeq,beq,lb,ub,x,options);
        xx=abs(xx);
        fxx=feval(fs,xx);
        np=np+1;
        nf=nf+1;
        diff=x-xx;
        flag=0;
        %control if xx is a solution
        if norm(diff) <= eps
            rx=diff;
            x=xx;
            if iter==0
                iter=1;
                alphavis(iter)=alpha;
            end; 
        else    
            fdiff=fx-fxx;
            num=fdiff'*fdiff;
            den=diff'*diff;
            if ((alpha^2)*(num/den))>(beta^2 )
                %adaptive stepleght 
                alpha=min((alpha/2),sqrt(den/(2*num)));
                flag=1;
            else  
                xold=x;
                f=-(x-alpha*fxx);
                options=optimset('LargeScale','off','display','off');
                [x]=quadprog(P,f,A,b,Aeq,beq,lb,ub,xx,options);  
                x=abs(x);
                fx=feval(fs,x);
                np=np+1;
                nf=nf+1;
                rx=xold-x;
                iter=iter+1;
                alphavis(iter)=alpha;   
            end;
        end;  
end; 
t=cputime-tempo
fprintf('solution=%g \n',x);
fprintf('number of projections=%g \n',np);          
fprintf('number of function evaluations=%g \n',nf);         
fprintf('number of iterations=%g \n',iter);
