N = 1e3;

tol = 4e-12;
j = 0;
for i = 1:N
  a = reshape (rand*[1:16], 4, 4);  % Rank 2 matrix
  b = pinv (a);
  if (norm(a*b*a - a)>tol)
     warning('a*b*a-a')
     j=j+1;
  end
  if (norm(b*a*b - b)>tol)
     warning('b*a*b-b')
     j=j+1;
  end
  if (norm ((b*a)'-b*a)>tol)
     warning('(b*a)\'-b*a')
     j=j+1;
  end
  if (norm ((a*b)'-a*b)>tol)
     waning('(a*b)\'-a*b')
     j=j+1;
  end
end
warning_percentage = j/(4*N)
