Hi,
I'm having trouble debugging a program that I'm writing to reduce a matrix into Reduced Row-Echelon form (RREF). I'm using recursion and Gaussian Elimination to accomplish this task.
The problem is that the last values of the Matrix 'M' in the recursive loop don't seem to carry outside of the recursion; that is, the operations that I perform on them do not seem to stick. This is easier to show with part of the debug output:
Getting out of the recursion, M is... 1.0000 0.8333 Check on m 1.0000 0.8333 At the end of the function, M is... 1.0000 0.8333 Getting out of the recursion, M is... 1 3 2 0 6 5
Full code of the Gaussian Elimination recursive program as below:
function [M,pivots] = echelon(M) %Reduce matrix to reduced row-echelon form
%ECHELON(M) Returns reduced row-echelon matrix M and pivots vector 'pivots'
%Check if input is numerical matrix. If it isn't, give error.
narginchk(1,1); if ~(ismatrix(M) && isnumeric(M)) error('Not a numeric array'); end %Check if the matrix is empty. If it is, return
pivots=[]; if isempty(M) disp('found empty'); pivots = []; return; end %Find size of Matrix, calculate tolerance level
[m,n] = size(M); tolerance = length(M)*norm(M)*eps; %Find the pivot variable - it is value x, at position p
[x,p] = max(abs(M(:,1))); if x>tolerance M = Exchange(M,1,p); M = Divide(M,1,M(1,1)); for i=2:m if M(i,1)==0 %do nothing
else Subtract(M,1,i,M(i,1)/M(1,1)) end end disp('Recurring now.'); [~,newpivots] = echelon(M(2:m,2:n)); disp('Getting out of the recursion, M is...'); disp(M); pivots(1) = 1; for z=1:length(newpivots) pivots(end+1) = newpivots(z)+1; end for z=2:length(pivots) currentpivot = pivots(z); M = Subtract(M,z,1,M(1,currentpivot)/M(z,currentpivot)) end fprintf('Check on m\n'); disp(M); else disp('the else has been called.'); [~,newpivots] = echelon(M(:,2:n)); for z=1:length(newpivots) pivots(end+1) = newpivots(z) + 1; end end disp('At the end of the function, M is...'); disp(M); end
If you're familiar with Gaussian row reduction, the subtract, exchange, and divide functions used in the program are all elementary row operations (subtracting row i2 by row i1 * a scalar).
Functions as follows:
function returnMatrix = Subtract(M,i1,i2,r)if ~(nargin == 4) error('# of arguments incorrect. Exiting.');endif ~((i1 >= 1) && (i1 <= size(M,1)) && (i2 >= 1) && (i2 <= size(M,1))) error('Row i not in bounds of matrix. Exiting');endif r==0 error('r is zero. Exiting.');endreturnMatrix = M;returnMatrix(i2,:) = returnMatrix(i2,:) - (r*returnMatrix(i1,:));endfunction returnMatrix = Exchange(M,i1,i2)if ~(nargin == 3) error('Incorrect arguments. Exiting.');endreturnMatrix = M;if (i1==round(i1)) && (i2==round(i2)) && (i1 >= 1) && (i2 >= 1) && (i1 <= size(M,1)) && (i2 <= size(M,1)) returnMatrix([i1 i2],:) = returnMatrix([i2 i1],:);else error('Incorrect arguments entered for row values. Exiting.');endendfunction returnMatrix = Divide(M,i,r)if ~(nargin == 3) error('# of arguments incorrect. Exiting.');endif ~((i >= 1) && (i <= size(M,1))) error('Row i not in bounds of matrix. Exiting');endif r==0 error('r value is zero. Exiting');endreturnMatrix = M;returnMatrix(i,:) = returnMatrix(i,:)/r;end
Sorry, this is a pretty big post, and I've been working on this small error for hours. I was wondering if someone can help me solve the problem I'm having?
Best Answer