function [W,IPIVOT,IFLAG] = factor_cd(W)

%*****************************************************************************80
%
%MATLAB version of SUBROUTINE FACTOR on pages 165-6 of conte/deboor
%**************   I N P U T   ************
%  W  array of size N by N containing the matrix  A  of order  N to be factored
%**************  O U T P U T  ************
%  W  array of size N by N containing the factorization of  P*A  for some
%     permutation matrix  P  specified by  IPIVOT.
%  IPIVOT  integer vector of length  N  indicated that row  IPIVOT(k)
%     was used to eliminate the k-th unknown, k=1:N.
%  IFLAG  an integer
%       =  1  if an even number of interchanges were carried out
%       = -1  if an off number of interchanges were carried out
%       =  0  if the upper triangular factor has one or more 
%               zero diagonal entries.
%     Thus, determinant(A) = IFLAG*W(1,1)*...*W(N,N)
%     IF IFLAG~=0, then the linear system  A*X = B is obtained as
%         X = SUBST(W,IPIVOT,B); 
%**************  M E T H O D  ************
%  The program follows Algorithm 4.2 using scaled partial pivoting.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%

% initialize IFLAG,  IPIVOT,  D

IFLAG = 1; 

[N,~] = size(W); if N==0, disp('input  W  is empty'), return, end

IPIVOT = 1:N; D = max(abs(W')); I = D == 0; 
if ~isempty(D(I)), IFLAG = 0; D(I) = 1; end 
if N == 1, return, end

for k=1:N-1
%            determine the pivot row, the row  ISTAR.
   [~,ISTAR] = max(abs(W(k:N,k))./D(k:N)); ISTAR = ISTAR + k-1;
   if W(ISTAR,k)==0, IFLAG = 0;
   else
      if ISTAR> k  % switch row k with row ISTAR to make k the pivot row
         IFLAG = -IFLAG;
         W([k,ISTAR],:) = W([ISTAR,k],:); 
         IPIVOT([k,ISTAR]) = IPIVOT([ISTAR,k]);
         D([k,ISTAR]) = D([ISTAR,k]);
      end
            % use row  k  to eliminate  X(k)  from rows k+1:N
      W(k+1:N,k) = W(k+1:N,k)/W(k,k);
      W(k+1:N,k+1:N) = W(k+1:N,k+1:N) - W(k+1:N,k)*W(k,k+1:N);
   end
end
if W(N,N) == 0; IFLAG = 0; end
return, end
