30 lines
523 B
Mathematica
30 lines
523 B
Mathematica
|
function x = solveLrPivot(A, P, b)
|
||
|
|
len = length(b)
|
||
|
|
z = zeros(len,1)
|
||
|
|
|
||
|
|
% solve Lz = b
|
||
|
|
for i = 1:len
|
||
|
|
z(i) = b(i)
|
||
|
|
|
||
|
|
for j = 1:(i-1)
|
||
|
|
z(i) = z(i) - z(j) * A(P(i), j)
|
||
|
|
end
|
||
|
|
|
||
|
|
% no division necessary as L is unipotent
|
||
|
|
end
|
||
|
|
|
||
|
|
x = zeros(len,1)
|
||
|
|
|
||
|
|
% solve Rx = z
|
||
|
|
for i = len:-1:1
|
||
|
|
x(i) = z(i)
|
||
|
|
|
||
|
|
for j = len:-1:(i+1)
|
||
|
|
x(i) = x(i) - x(j) * A(P(i), j)
|
||
|
|
end
|
||
|
|
|
||
|
|
x(i) = x(i) / A(P(i), i)
|
||
|
|
end
|
||
|
|
|
||
|
|
x_ret = x
|
||
|
|
end
|