function [dx,dy,ds] = newtonsolve(NS,r,p,q)
%NEWTONSOLVE    Solve linear system with factorized matrix.
%
%  [dx,dy,ds] = newtonsolve(NS,r,p,q)

  m = size(r,1);
  n = size(p,1);

  % ********************************************************************
  % Solve KKT system with LU factors.
  % ********************************************************************

  if (NS.method == 1)
    rhs = [p-q./NS.x; r];
    warn_stat = warning;
    warning('off','all');
    lhs = NS.U\(NS.L\rhs(NS.pp));
    warning(warn_stat);
    lhs(NS.qq) = lhs;
    dx = lhs(1:n);
    dy = lhs(n+1:n+m);
    ds = (q-NS.s.*dx)./NS.x;
  end


  % ********************************************************************
  % Solve KKT system with LDL' factors.
  % ********************************************************************

  if (NS.method == 2)
    rhs = [p-q./NS.x; r];
    warn_stat = warning;
    warning('off','all');
    lhs = NS.L'\(NS.D\(NS.L\rhs(NS.pp)));
    warning(warn_stat);
    lhs(NS.pp) = lhs;
    dx = lhs(1:n);
    dy = lhs(n+1:n+m);
    ds = (q-NS.s.*dx)./NS.x;
  end
end
