% This computes the search direction for IPM_DK.
%
% Jan 25, 2012

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

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

if r < m
  
  d0 = w.^2;
  tmp = full(2*skron(blk,invX,KXinvWKX) - skron(blk,KXinvWKX,KXinvWKX));
  tmp = (tmp + tmp')/2;
  [V,D] = svd(tmp);
  V = V(:,1:r);
  D = D(1:r,1:r);
  
  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,r+1,r+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(2*skron(blk,invX,KXinvWKX) - skron(blk,KXinvWKX,KXinvWKX));
  tmp = (tmp + tmp')/2;
  
  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