This question contains edits
I'm currently trying to obtain a PA=LU factorisation with Octave. Here is my algorythm:
The main script:
clear all; close all; clc;
%{
We have a square matrix of order 5.
We can show that L4*P4*L3*P3*L2*P2*L1*P1*A = U
With L'4 = L4
L'3 = P2*L3*P4
L'2 = P4*P3*L2*P3*P4
L'1 = P4*P3*P2*L1*P2*P3*P4
and L = \prod L'i^-1
P = P4*P3*P2*P1
Then, for n = 5, we find Pi and Li with:
Pi = eye(n);
Pi([i, n], i) = Pi([n, i], i);
Li = eye(5);
Li(i 1:n, i) = -A(i 1:n, i)/A(i, i); % We divide by the pivot
We will write a function "pli(A, i, n, LorP)" for that.
%}
A = [2 6 5 1 7; ...
6 9 2 2 3; ...
7 9 0 6 5; ...
1 1 6 6 7; ...
2 7 7 3 9]; % the given matrix
Asv = A; % because we will change the value of A
P1 = pli(A, 1, 5, "P")
A = P1*A;
L1 = pli(A, 1, 5, "L")
A = L1*A;
P2 = pli(A, 2, 5, "P")
A = P2*A;
L2 = pli(A, 2, 5, "L")
A = L2*A;
P3 = pli(A, 3, 5, "P")
A = P3*A;
L3 = pli(A, 3, 5, "L")
A = L3*A;
P4 = pli(A, 4, 5, "P")
A = P4*A;
L4 = pli(A, 4, 5, "L")
A = L4*A;
P = P4*P3*P2*P1
Lp4 = L4;
Lp3 = P2*L3*P4;
Lp2 = P4*P3*L2*P3*P4;
Lp1 = P4*P3*P2*L1*P2*P3*P4;
L = eye(5);
L(1:5, 1) = Lp1(1:5, 1);
L(2:5, 2) = Lp2(2:5, 2);
L(3:5, 3) = Lp3(3:5, 3);
L(4:5, 4) = Lp4(4:5, 4);
L
U = A
% Verification:
[l u p] = lu(Asv);
norm(p-P)
norm(l-L)
norm(u-U)
The function pli(A, i, n, LorP):
function Mi = pli(A, i, n, LorP)
% order n
if LorP == "P"
Pi = eye(n);
Pi([i, n], i) = Pi([n, i], i); % matrix of elementary permutations
Mi = Pi;
elseif LorP == "L"
Li = eye(5);
Li(i 1:n, i) = -A(i 1:n, i)/A(i, i); % Devide by the pivot
Mi = Li;
else
error("Wrong entry for 'LorP'")
endif
endfunction
When I execute my code, I obtain strange infinite values, and don't understand why, and how to solve it:
P1 =
0 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
1 0 0 0 1
L1 =
1 0 0 0 0
-Inf 1 0 0 0
-Inf 0 1 0 0
-Inf 0 0 1 0
-Inf 0 0 0 1
P2 =
1 0 0 0 0
0 0 0 0 0
0 0 1 0 0
0 0 0 1 0
0 1 0 0 1
L2 =
1 0 0 0 0
0 1 0 0 0
0 NaN 1 0 0
0 NaN 0 1 0
0 NaN 0 0 1
P3 =
1 0 0 0 0
0 1 0 0 0
0 0 0 0 0
0 0 0 1 0
0 0 1 0 1
L3 =
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 NaN 1 0
0 0 NaN 0 1
P4 =
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 0 0
0 0 0 1 1
L4 =
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 NaN 1
P =
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
1 1 1 1 1
L =
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
U =
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
ans = 2.2361
DLASCL: parameter number 4 is invalid
error: Fortran procedure terminated by call to XERBLA
error: called from
q1a at line 69 column 1
Thank you in advance for your help.
EDIT (related to the answer of @BillBokeey)
After changing pli.m into:
function Mi = pli(A, i, n, LorP)
% order n
if LorP == "P"
Pi = eye(n);
if(A(i, i)==0)
idx_Pvt = find(diag(A)~=0,1); % Find index of first nonzero pivot
Pi([i, idx_Pvt], [i,idx_Pvt]) = [0 1;1 0];
end
Mi = Pi;
elseif LorP == "L"
Li = eye(5);
Li(i 1:n, i) = -A(i 1:n, i)/A(i, i); % Devide by the pivot
Mi = Li;
else
error("Wrong entry for 'LorP'")
endif
endfunction
the answer becomes:
P1 =
Diagonal Matrix
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
L1 =
1.000000000000000 0 0 0 0
-3.000000000000000 1.000000000000000 0 0 0
-3.500000000000000 0 1.000000000000000 0 0
-0.500000000000000 0 0 1.000000000000000 0
-1.000000000000000 0 0 0 1.000000000000000
P2 =
Diagonal Matrix
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
L2 =
1.000000000000000 0 0 0 0
0 1.000000000000000 0 0 0
0 -1.333333333333333 1.000000000000000 0 0
0 -0.222222222222222 0 1.000000000000000 0
0 0.111111111111111 0 0 1.000000000000000
P3 =
Diagonal Matrix
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
L3 =
1.000000000000000e 00 0 0 0 0
0 1.000000000000000e 00 0 0 0
0 0 1.000000000000000e 00 0 0
0 0 3.833333333333307e 01 1.000000000000000e 00 0
0 0 3.333333333333310e 00 0 1.000000000000000e 00
P4 =
Diagonal Matrix
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
L4 =
1.000000000000000 0 0 0 0
0 1.000000000000000 0 0 0
0 0 1.000000000000000 0 0
0 0 0 1.000000000000000 0
0 0 0 -0.096069868995633 1.000000000000000
P =
Diagonal Matrix
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
L =
1.000000000000000e 00 0 0 0 0
-3.000000000000000e 00 1.000000000000000e 00 0 0 0
-3.500000000000000e 00 -1.333333333333333e 00 1.000000000000000e 00 0 0
-5.000000000000000e-01 -2.222222222222222e-01 3.833333333333307e 01 1.000000000000000e 00 0
-1.000000000000000e 00 1.111111111111111e-01 3.333333333333310e 00 -9.606986899563325e-02 1.000000000000000e 00
U =
2.000000000000000e 00 6.000000000000000e 00 5.000000000000000e 00 1.000000000000000e 00 7.000000000000000e 00
0 -9.000000000000000e 00 -1.300000000000000e 01 -1.000000000000000e 00 -1.800000000000000e 01
0 0 -1.666666666666679e-01 3.833333333333333e 00 4.500000000000000e 00
0 0 0 1.526666666666656e 02 1.799999999999988e 02
0 0 0 1.776356839400250e-15 -2.292576419213974e 00
ans = 1.902113032590308
ans = 38.52133878794479
ans = 241.8726657311275
CodePudding user response：
The main problem in your algorithm is that your matrix P is not a permutation matrix (e. g. for P1, You have no 1s on the first row and two 1s on the last row).
Furthermore, you don't necessarily need to permute your matrix A at every step, just when the pivot is nil. This is the reason why the permutation matrix is here in the first place
Both these problem can be solved by modifying your function pli :
function Mi = pli(A, i, n, LorP)
% order n
if strcmp(LorP,"P")
Pi = eye(n);
if(A(i, i)==0)
idx_Pvt = find(diag(A)~=0,1); % Find index of first nonzero pivot
Pi([i, idx_Pvt], [i,idx_Pvt]) = [0 1;1 0];
end
Mi = Pi;
elseif strcmp(LorP,"L")
Li = eye(5);
Li(i 1:n, i) = -A(i 1:n, i)/A(i, i); % Devide by the pivot
Mi = Li;
else
error("Wrong entry for 'LorP'")
endif
endfunction
Like I said in one of my comment, you also need to modify your main in order to build the L ans U matrices correctly. Here is a (non-optimised) MATLAB example of how to do this, which should work in Octave after a small amount of modifications:
A = [2 6 5 1 7; ...
6 9 2 2 3; ...
7 9 0 6 5; ...
1 1 6 6 7; ...
2 7 7 3 9]; % the given matrix
[sz1,~] = size(A);
A_ite = A;
L_out = eye(sz1);
L_Save = zeros(sz1,sz1,sz1);
for ii=1:sz1
% Calculate Pi
Pi = pli(A_ite, ii, sz1, 'P');
A_ite = Pi*A_ite;
Li = pli(A_ite,ii,sz1,'L');
A_ite = Li*A_ite;
L_Save(:,:,ii) = (Li-eye(sz1));
L_out = Li*Pi*L_out;
end
U = A_ite;
L = -sum(L_Save,3) eye(sz1);
norm(A-L*U)
You can also find an efficient version of the algorithm at the end of this page