Inversion of a matrix by pivoting method

pivinv.m

function [A_inv, d] = pivinv(A)
% PIVINV inversion of a matrix by pivoting method
%
% Input: A, a square matrix
% Outputs: A_inv, the inverse of A; d, the determinant of A
%
% Reference: Enrique Castillo and Francisco Jubete, "The Gamma-algorithm
% and some applications," International Journal of Mathematical Education
% in Science and Technology 35, no. 3 (2004): 369-389.
%
% Coded by steven2358 in 2006.

if (size(A,1) ~= size(A,2)),
    error('Matrix must be square');
end

N = size(A,1);		% matrix order
V = cell(N,1);		% intermediate matrices
p = cell(N,1);		% scalar products
pivs = zeros(N,1);	% pivots

V_prev = eye(N);
V_new = zeros(N);
for i=1:N,
    A_row = A(i,:);
    
    p{i} = A_row*V_prev; % scalar products (pivots)
    
    pivind = find(p{i}(i:end)); % index of first pivot, different from zero
    if isempty(pivind),
        error('Singular matrix, range %d',i);
    else
        pivind = pivind(1) + i - 1;
    end
    pivs(i) = p{i}(pivind);
    
    pivcol = V_prev(:,pivind)/p{i}(pivind); % pivoting transformation
    for j=1:N,
        V_new(:,j) = V_prev(:,j)-p{i}(j)*pivcol;
    end
    V_new(:,pivind) = pivcol;
    V{i} = V_new;
    V_prev = V_new;
end

d = sum(pivs);
A_inv = V{N};