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