% This computes the search direction for IPM_p.
%
% Sep 16, 2012

function dir = search_direction(w,MA,blk,ME,p,mu,g,V,diagD,dir_in)

n = blk{1,2};
nn = n*(n+1)/2;
m = size(w,1);

dg = zeros(nn,1);
count = 1;
for i=1:n
  for j=1:i
    if diagD(i) ~= diagD(j)
      dg(count) = p*(diagD(i)^(p-1)-diagD(j)^(p-1))/(diagD(i)-diagD(j));
    else
      dg(count) = p*(p-1)*diagD(i)^(p-2);
    end
    count = count + 1;
  end
end

kronV = full(skron(blk,V,V));

if nn < m
  
  d0 = w.^2;
  tmp = kronV*spdiags(dg,0,nn,nn)*kronV';  
  [V,D] = svd(tmp);
  
  V = [V' zeros(size(V,2),1); zeros(1,size(V,1)) 1];
  d = [mu./diag(D); (1 - sum(w))^2];   
  R = spdiags(d,0,nn+1,nn+1) + V*(ME*spdiags(d0,0,m,m)*ME')*V';
  
  [U,D] = svd(R);
  d2 = diag(D).^(-1);
  gd = g.*d0;
  dir = 1/mu*(-gd + (ME'*(V'*(U*(d2.*(U'*(V*(ME*gd))))))).*d0);  
  
else
  
  tmp = kronV*spdiags(dg,0,nn,nn)*kronV';
  
  tmp = [tmp zeros(size(tmp,1),1);zeros(1,size(tmp,1)) mu/(1 - sum(w))^2];
  d0 = mu./(w.^2);
  
  Hhandle = @(x) Hessfun(x,ME,tmp,d0);
  
  [dir, dummy] = pcg(Hhandle,-g,1e-10,100,[],[],dir_in);  
  
end

%%% Nested function

function y = Hessfun(x,ME,tmp,d0)

  y = ME'*(tmp*(ME*x)) + d0.*x;
    
return