I'm trying to create a program that takes a square (n-by-n) matrix as input, and if it is invertible, will LU decompose the matrix using Gaussian Elimination.
Here is my problem: in class we learned that it is better to change rows so that your pivot is always the largest number (in absolute value) in its column. For example, if the matrix was
A = [1,2;3,4] then switching rows it is
[3,4;1,2] and then we can proceed with the Gaussian elimination.
My code works properly for matrices that don't require row changes, but for ones that do, it does not. This is my code:
function newgauss(A) [rows,columns]=size(A); P=eye(rows,columns); %P is permutation matrix if(det(A)==0) %% determinante is 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 for i=pivot+1:rows if(abs(U(i,pivot))>max) max=abs(U(i,pivot)); maxi=i; end end %%if needed then 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 the column pivot using gauss elimination for 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 we are switching rows, we are looking to decompose
P (permutation matrix) times
A, and not the original
A that we had as input.
Explanation of the code:
- First I check if the matrix is invertible, if it isn't, stop. If it is, pivot is (1,1)
- I find the largest number in column 1, and switch rows
- Grade column 1 using Gaussian elimination, turning all but the spot (1,1) to zero
- Pivot is now (2,2), find largest number in column 2... Rinse, repeat