% This computes the search direction for IPM_D.
%
% Jan 12, 2012

function dir = search_direction(X,blk,w,MA,ME,mu,g,invX,dir_in)

m = size(MA,2);
n = size(X,1);
nn = n*(n + 1)/2;

if nn < m
  
  d0 = w.^2;
  tmp = full(skron(blk,invX,invX));
  [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 = full(skron(blk,invX,invX));
  
  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