LU decomposing a square matrix matlab gauss elimination

Question

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:

1. First I check if the matrix is invertible, if it isn't, stop. If it is, pivot is (1,1)
2. I find the largest number in column 1, and switch rows
3. Grade column 1 using Gaussian elimination, turning all but the spot (1,1) to zero
4. Pivot is now (2,2), find largest number in column 2... Rinse, repeat

Answer 1

Your code seems to work fine from what I can tell, at least for the basic examples A=[1,2;3,4] or A=[3,4;1,2]. Change your function definition to:

function [L,U,P] = newgauss(A)

so you can output your calculated values (much better than using disp, but this shows the correct results too). Then you'll see that P*A = L*U. Maybe you were expecting L*U to equal A directly? You can also confirm that you are correct via Matlab's lu function:

[L,U,P] = lu(A);
L*U
P*A

Permutation matrices are orthogonal matrices, so P⁻¹ = P^T. If you want to get back A in your code, you can do:

P'*L*U

Similarly, using Matlab's lu with the permutation matrix output, you can do:

[L,U,P] = lu(A);
P'*L*U

(You should also use error or warning rather than how you're using disp in checking the determinant, but they probably don't teach that.)