algorithm - LU decomposing a square matrix matlab gauss elimination -

i'm trying create program takes square (n-by-n) matrix input, , if invertible, lu decompose matrix using gaussian elimination.

here problem: in class learned better change rows pivot largest number (in absolute value) in column. example, if matrix a = [1,2;3,4] switching rows [3,4;1,2] , can proceed gaussian elimination.

my code works matrices don't require row changes, ones do, not. code:

function newgauss(a)     [rows,columns]=size(a);     p=eye(rows,columns); %p permutation matrix     if(det(a)==0) %% determinante 0 means no single solution         disp('no solutions or infinite number of solutions')         return;     end     u=a;     l=eye(rows,columns);     pivot=1;     while(pivot<rows)         max=abs(u(pivot,pivot));         maxi=0;%%find maximum abs value in column pivot         i=pivot+1:rows             if(abs(u(i,pivot))>max)                 max=abs(u(i,pivot));                 maxi=i;             end         end %%if needed switch         if(maxi~=0)             temp=u(pivot,:);             u(pivot,:)=u(maxi,:);             u(maxi,:)=temp;             temp=p(pivot,:);             p(pivot,:)=p(maxi,:);             p(maxi,:)=temp;         end %%grade column pivot using gauss elimination         i=pivot+1:rows             num=u(i,pivot)/u(pivot,pivot);             u(i,:)=u(i,:)-num*u(pivot,:);             l(i,pivot)=num;         end         pivot=pivot+1;     end     disp('pa is:');     disp(p*a);     disp('lu is:');     disp(l*u); end 

clarification: since switching rows, looking decompose p (permutation matrix) times a, , not original a had input.

explanation of code:

1. first check if matrix invertible, if isn't, stop. if is, pivot (1,1)
2. i find largest number in column 1, , switch rows
3. grade column 1 using gaussian elimination, turning spot (1,1) zero
4. pivot (2,2), find largest number in column 2... rinse, repeat

your code seems work fine can tell, @ least basic examples a=[1,2;3,4] or a=[3,4;1,2]. change function definition to:

function [l,u,p] = newgauss(a) 

so can output calculated values (much better using disp, shows correct results too). you'll see p*a = l*u. maybe expecting l*u equal a directly? can confirm correct via matlab's lu function:

[l,u,p] = lu(a); l*u p*a 

permutation matrices orthogonal matrices, p⁻¹ = p^t. if want a in code, can do:

p'*l*u 

similarly, using matlab's lu permutation matrix output, can do:

[l,u,p] = lu(a); p'*l*u 

(you should use error or warning rather how you're using disp in checking determinant, don't teach that.)