Variable not advancing as told to

I'm attempting to create MATLAB code that transforms any matrix into its row reduced echelon form and logs the steps to LaTeX.

This is what I've done so far (just the steps to row echelon form):

function[lA] = rrefLaTeX(A)
%set initial parameters (j= index of leftmost non-zero column)
disp(A)
[m,n]=size(A);
i=1
j=find(any(A),1)
%if A(i,j)=0, switch row 1 with row k, where k is the index of the first row with a non-zero element in column j
while i<m
  if ((A(i,j) == 0))
    k = find(A(:,j),1);
    A([j,k],:)=A([k,j],:);
    disp('\overleftrightarrow{{l_j}\leftrightarrow{l_k}}')
    disp(A)
  end
  %make all elements below the pivot of row i equal 0
  for p = (i+1): m
    if logical(A(p,j)==0)==0
      A(p,:)=A(p,:)-(A(p,j)/A(i,j))*A(i,:);
      disp('\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}')
      disp(A)
    end
  end
  i=i+1
  j=find(any(A([i:end],:)),1)
end
end

This is what I obtain when I run it (with a random matrix A):

if true
  >> rrefLaTeX(A)
        3         -2         -5         -1
       -1          1          5          4
        4         -2          5         -2
       -5          0         -2          3
i = 1
j = 1
\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}
        3         -2         -5         -1
        0        1/3       10/3       11/3
        4         -2          5         -2
       -5          0         -2          3
\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}
        3         -2         -5         -1
        0        1/3       10/3       11/3
        0        2/3       35/3       -2/3
       -5          0         -2          3
\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}
        3         -2         -5         -1
        0        1/3       10/3       11/3
        0        2/3       35/3       -2/3
        0      -10/3      -31/3        4/3
i = 2
j = 2
\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}
        3         -2         -5         -1
        0        1/3       10/3       11/3
        0          0          5         -8
        0      -10/3      -31/3        4/3
\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}
        3         -2         -5         -1
        0        1/3       10/3       11/3
        0          0          5         -8
        0          0         23         38
i = 3
j = 2
\overleftrightarrow{{l_j}\leftrightarrow{l_k}}
        0        1/3       10/3       11/3
        3         -2         -5         -1
        0          0          5         -8
        0          0         23         38
warning: division by zero
warning: called from
  rrefLaTeX at line 20 column 11
\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}
        0        1/3       10/3       11/3
        3         -2         -5         -1
        0          0          5         -8
      0/0        0/0        1/0        1/0
i = 4
j = 1
end

The problem is that, when i changes to 3, j should change to 3 as well but it's not happening, making it all fail afterwards. Can someone explain what is going on?

Answer 1

[dpb]

"...when i changes to 3, j should change to 3 as well but it's not happening"

Looks like perfect place to use debugger to see where your logic error lies...set a breakpoint and trace operation.

The only thing that modifies j is

i=i+1
j=find(any(A([i:end],:)),1)

so the logic flaw has to be in this result; the condition you think is so isn't--the first element of the next row isn't zero; is it possible it isn't quite zero but small enough that disp shows it as such?

Since you make modifications to A by

A(p,:)=A(p,:)-(A(p,j)/A(i,j))*A(i,:);

it's quite possible I'd think that floating point roundoff is biting you (in fact, I'd say it's probable).

[Rodrigo Pereira]

How do I fix it then? How do I make it exactly zero?

[dpb]

Did you confirm that's the problem for certain?

You can only make it identically zero by testing and setting it for some tolerance; generally use a tolerance for the test and live with the result. Welcome to the wonderful and sometimes apparently wacky world of floating point--- :)

<Why is 0.3-0.2-0.1 not equal to zero?> is apropos to the question/problem.

if abs(x)<3*eps(x), x=0; end  % should give the idea

I thought there was a "fuzzy" comparison in later revisions, but a search in the doc didn't find what it is named...maybe one of the other regulars who recalls will come along and point it out.

[Rodrigo Pereira]

I modified my code. What I think I've done is: I turned the "almost zeroes" into zeroes before switching rows and before and after the row transformation.(After the transformation, I made sure that all zeroes of the row were exact). Yet, it still does not work, when i=3, j=2 and not 3 as it should.

if true
 function[lA] = rrefLaTeX(A)
%set initial parameters (j= index of leftmost non-zero column)
disp(A)
[m,n]=size(A);
i=1
j=find(any(A),1)

while i<m
%if A(i,j)=0, switch row 1 with row k, where k is the index of the first row with a non-zero element in column j
if (abs(A(i,j)) < eps)
  A(i,j)=0
  k = find(any(A([i+1:end],j),1))+i;
  disp(['k=',int2str(k)])
  A([i,k],:)=A([k,i],:);
  disp('\overleftrightarrow{{l_j}\leftrightarrow{l_k}}')
  disp(A)
end

%make all elements below the pivot of row i equal 0
for p = (i+1): m
  if (abs(A(p,j)) < eps*(A(p,j)))
  A(p,j)=0
  else
  A(p,:)=A(p,:)-(A(p,j)/A(i,j))*A(i,:);
  for a=1:n
  if (abs(A(p,a)) < eps*(A(p,a)))
  A(p,a)=0
  end
  end
  disp('\overleftrightarrow{{l_p}\gets{l_p-((A(p,j)/A(i,j))l_i}}')
  disp(A)
  disp(['p=',int2str(p)])
  end
end

i=i+1
j=find(any(A([i:end],:)),1)
end

end

end

[dpb]

Use the debugger to see where the logic error is...it's often the case when you actually step thru the code there's a real "AHA!" moment as one realizes what happened...

Did you confirm the issue about floating point rounding was really an issue; you never mentioned on that point...