%gradient projection method for max { g(x,u) | || u || <= rho }.
%Algorithm 1.3.16 in Polak-book
%set up for rho = 1 and delta = 1
%function y = psi(x)
%requires initial vectors x
%requires functions g(x,u), gradgu(x,u)
%returns matrix H containing iterates and values.
%Copyright Elijah Polak 03-14-98
function y = psi(x)
rho = 1;
delta = 1;
u = [1 1 1 1]';
u = rho*u/norm(u);
k=0;
u = u(:);
x = x(:);
F = -g(x,u);
G = -gradgu(x,u);
unew = rho*(delta*u - G)/norm(delta*u - G);
h = unew - u;
theta = G'*h + 0.5*delta*h'*h;
while theta < -1e-9
  h = rho*(delta*u - G)/norm(delta*u - G) - u;
  h = h(:);
    if -g(x,u + 0.8^k*h) - F >  0.8^k*0.8*theta
     while -g(x,u + 0.8^k*h) - F >  0.8^k*0.8*theta
       k=k+1;
         if k > 50, break
         end %if
     end %while
    else
     while -g(x,u + 0.8^k*h) - F <= - 0.8^k*0.8*h'*h
        k=k-1;
                if k < 0, break
                end %if
        end %while
        k=k+1;
  end %if
  u = u + 0.8^k*h;
  F = -g(x,u);
  G = -gradgu(x,u);
  h = rho*(delta*u - G)/norm(delta*u - G) - u;
  theta = G'*h + 0.5*delta*h'*h;
end %while
y = g(x,u);