Could MATLAB please provide rref() without the rational conversion?

I am having the same problem as described in:

https://au.mathworks.com/matlabcentral/answers/102622-why-does-rref-give-significantly-different-answers-than-for-solutions-of-some-systems-in-matlab

The answer given is an attached MATLAB file, "rrefNoRational.m" which is a copy of "rref.m" without the code that performs the rational conversion that is the cause of the problem. However, I cannot find the attached file. Could the file rrefNoRational.m please be provided?

Thank you for your time.

Answer 1

[John D'Errico]

Just look at the code.

type rref

function [A, jb] = rref(A, tol)
%RREF   Reduced row echelon form.
%   R = RREF(A) produces the reduced row echelon form of A.
%
%   [R,jb] = RREF(A) also returns a vector, jb, so that:
%       r = length(jb) is this algorithm's idea of the rank of A,
%       x(jb) are the bound variables in a linear system, Ax = b,
%       A(:,jb) is a basis for the range of A,
%       R(1:r,jb) is the r-by-r identity matrix.
%
%   [R,jb] = RREF(A,TOL) uses the given tolerance in the rank tests.
%
%   Roundoff errors may cause this algorithm to compute a different
%   value for the rank than RANK, ORTH and NULL.
%
%   Class support for input A:
%      float: double, single
%
%   See also RANK, ORTH, NULL, QR, SVD.

%   Copyright 1984-2017 The MathWorks, Inc. 

[m,n] = size(A);

% Does it appear that elements of A are ratios of small integers?
[num, den] = rat(A);
rats = isequal(A, num./den);

% Compute the default tolerance if none was provided.
if (nargin < 2)
    tol = max(m,n)*eps(class(A))*norm(A,inf);
end

% Loop over the entire matrix.
i = 1;
j = 1;
jb = zeros(1,0);
while i <= m && j <= n
   % Find value and index of largest element in the remainder of column j.
   [p, k] = max(abs(A(i:m,j)));
   k = k+i-1;
   if p <= tol
      % The column is negligible, zero it out.
      A(i:m,j) = 0;
      j = j + 1;
   else
      % Remember column index
      jb = [jb j]; %#ok<AGROW>
      % Swap i-th and k-th rows.
      A([i k],j:n) = A([k i],j:n);
      % Divide the pivot row by the pivot element.
      A(i,j:n) = A(i,j:n)./A(i,j);
      % Subtract multiples of the pivot row from all the other rows.
      for k = [1:i-1 i+1:m]
         A(k,j:n) = A(k,j:n) - A(k,j).*A(i,j:n);
      end
      i = i + 1;
      j = j + 1;
   end
end

% Return "rational" numbers if appropriate.
if rats
    [num, den] = rat(A);
    A = num./den;
end

Do you see the very last thing it does is do the conversion? Just write your own version that does not have that last if clause in it, copied, except for those last 4 lines.