mardi 22 octobre 2019

How do I implement partial pivoting correctly within my PA = LU factorization algorithm?

I am working through some practice exercises for MatLab practice and one of the problems gives me pseudocode to follow:

STEP 1: Set P = (1:n)', set M = A

STEP 2: For i = 1,2,...,n-1 do STEPS 3-6 

STEP 3: Find k so that abs(M(k,i)) >= abs(M(l,i)) for l = i,i+1,...,n

STEP 4: Interchange rows i and rows k of M and entries i and k of P

STEP 5: If M(i,i) = 0
       OUTPUT("Matrix is rank deficient")

STEP 6: For j = i+1,i+2,...,n do STEPS 7-9

STEP 7: Set  a = M(j,i)/M(i,i)

STEP 8: Set  M(j,i:end) = M(j,i:end) - a*M(i,i:end)

STEP 9: Set  M(j,i) = a

STEP 10: OUTPUT([P,M]); STOP.

At the end of STEP 6, for each i, I have a matrix [P(i), A(i)] and I would need to construct the matrix [[P(1), A(1)], [P(2), A(2)], ... ,[P(n-1), A(n-1)]] (hence the use of vertcat) which will be a 30x7 matrix. However, I have trouble implementing partial pivoting within my algorithm and I don't end up with a 30x7 matrix. Here's what I have so far:

r = RandStream('mt19937ar','Seed',1234);
A = r.randn(6,6);
n = 6;
P = (1:n)';
M = A;
t1_out = [];

for ii = 1:n-1
    [~,k] = max(abs(M(ii:n,ii)));
    temp1 = M(ii,:);
    A(ii,:) = M(k + (ii-1),:);
    M(k + (ii-1),:) = temp1;
    temp2 = P(ii,1);
    P(ii,1) = P(k + (ii-1),1);
    P(k + (ii-1),1) = temp2;

    if M(ii,ii) == 0
        error('Matrix is rank deficient');
    end

    for j = ii+1:n
        a = M(j,ii)/M(ii,ii);
        M(j,ii:end) = M(j,ii:end) - a*M(ii,ii:end);
        M(j,ii) = a;
    end

    t1_out = vertcat(t1_out,[P,M]);
end

I try to find k by looking for the maximum value within the columns. I know the issue is somewhere within my partial pivoting code, but I don't what it is that causes it to output incorrectly. Any help or feedback is appreciated.

Aucun commentaire:

Enregistrer un commentaire