format compact
n=6;
%
% Create n x n linear system  A*x = b:
%
A=floor(10*rand(n))
b=-(1:n)'
%
% Random permutation:
%
p = randperm(n)
%
% Inverse permutation:
%
invp = zeros(1,n);
invp(p)=1:n
%
% Permutation matrices (for demo purposes only, we never form these
% matrices in practise):
%
P = eye(n);
P=P(p,:)
Pinv=P'
%
% Solve permuted linear system P*A*inv(P)*y=P*b:
%
y = A(p,p)\b(p);
%
% Get the solution to the original problem  x=inv(P)*y:
%
x = y(invp);
%
% Check that we really solved the original A*x=b:
%
fprintf(1,'The norm of difference between solutions to ');
fprintf(1,'P*A*inv(P)*P*x=P*b and A*x=b is %e \n',norm(x-A\b));
