Hi Ye,
The missing logic in your code is the part that identifies the pivot element needed for partial pivoting. Specifically, the steps to determine the pivot index 'k' are absent. Here's how you can implement it:
- Locate the Pivot: Before the if j~=k statement, you need to find the row index 'k' of the largest absolute value in the current column from row 'j' to 'n'.
- Implement Partial Pivoting: Use the pivot index k to swap rows in both 'U' and 'L' matrices if necessary.
Below is the complete implementation for LUP decomposition:
function [L, U, P] = mylu(A)
[m, n] = size(A);
if m ~= n
fprintf('The input matrix is not square!');
L = []; U = []; P = []; return;
end
L = eye(n, n);
U = A;
p = 1:n;
for j = 1:n-1
% Locate the pivot
% find row index k with the largest absolute
% value in the jth column
[~, k] = max(abs(U(j:n, j)));
k = k + j - 1; % Adjust index
if j ~= k
U([j k], j:n) = U([k j], j:n);
L([j k], 1:j-1) = L([k j], 1:j-1);
p([j k]) = p([k j]);
end
for i = j+1:n
L(i, j) = U(i, j) / U(j, j);
U(i, j) = 0;
U(i, j+1:n) = U(i, j+1:n) - L(i, j) * U(j, j+1:n);
end
end
P = eye(n);
P = P(p, :);
end
