function NS = newtonsolver(A,x,s)
%NEWTONSOLVER    Prepare linear solver for Newton system.
%
%  NS = newtonsolver(A,x,s) returns a MATLAB struct that holds the
%  factorization of the Newton matrix for solving the step equations in
%  the IPM.

  [m n] = size(A);
  NS = struct();

  % Factorization Method
  % --------------------
  % 1: augmented system, LU factorization
  % 2: augmented system, LDL' factorization (MATLAB only)
  % 3: normal equations, Cholesky factorization (not implemented)

  NS.method = 1;

  % ********************************************************************
  % LU factorization of KKT matrix
  % ==============================
  %
  % Reduce Newton system
  %
  %   [ A 0  0 ] (dx)   (r)
  %   [ 0 A' I ] (dy) = (p)
  %   [ S 0  X ] (ds)   (q)
  %
  % to the KKT system
  %
  %   [ -D^(-2) A' ] (dx) = (p-X\q)
  %   [  A      0  ] (dy)   (r),
  %
  % where D is diagonal with D(i,i) = sqrt(x(i)/s(i)). Denote the matrix
  % by K and calculate its LU factorization
  %
  %   K(pp,qq) = L*U.
  %
  % This factorization is available in both MATLAB and OCTAVE and uses
  % UMFPACK when K is sparse.
  % ********************************************************************

  if (NS.method == 1)
    NS.x = x;
    NS.s = s;
    K = [-diag(s./x), A'; A, sparse(m,m)];
    [NS.L,NS.U,NS.pp,NS.qq] = lu(K,'vector');
  end

  % ********************************************************************
  % LDL' factorization of KKT matrix
  % ================================
  %
  % Perform the same reduction as above but calculate its symmetric
  % indefinite factorization
  %
  %   K(pp,pp) = L*D*L'.
  %
  % MATLAB uses MA57 when K is sparse, which is not available in OCTAVE. 
  % ********************************************************************

  if (NS.method == 2)
    NS.x = x;
    NS.s = s;
    K = [-diag(s./x), A'; A, sparse(m,m)];
    [NS.L,NS.D,NS.pp] = ldl(K,'vector');
  end

end
