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

function dir = search_direction(X,r,w,MA,MB,ME,p,K,VtKtinvX,VDV,invX,mu,g,diagD,dir_in,blk)

k = size(K,2);
m = size(w,1);
kk = k^2;

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

kronVKX = kron(VtKtinvX,VtKtinvX);
pinvX = -invX;
kVMB = kronVKX*MB;

if r < m

  d0 = w.^2;
  tmp = 2*skron(blk,pinvX,VDV) + kVMB'*spdiags(dg,0,kk,kk)*kVMB;
  [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 = 2*skron(blk,pinvX,VDV) + MB'*(kronVKX'*spdiags(dg,0,kk,kk)*kronVKX)*MB;
  
  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