78 lines
2.5 KiB
Ada
78 lines
2.5 KiB
Ada
package body Decomposition is
|
|
|
|
procedure Swap_Rows (M : in out Matrix.Real_Matrix; From, To : Natural) is
|
|
Temporary : Matrix.Real;
|
|
begin
|
|
if From = To then
|
|
return;
|
|
end if;
|
|
for I in M'Range (2) loop
|
|
Temporary := M (M'First (1) + From, I);
|
|
M (M'First (1) + From, I) := M (M'First (1) + To, I);
|
|
M (M'First (1) + To, I) := Temporary;
|
|
end loop;
|
|
end Swap_Rows;
|
|
|
|
function Pivoting_Matrix
|
|
(M : Matrix.Real_Matrix)
|
|
return Matrix.Real_Matrix
|
|
is
|
|
use type Matrix.Real;
|
|
Order : constant Positive := M'Length (1);
|
|
Result : Matrix.Real_Matrix := Matrix.Unit_Matrix (Order);
|
|
Max : Matrix.Real;
|
|
Row : Natural;
|
|
begin
|
|
for J in 0 .. Order - 1 loop
|
|
Max := M (M'First (1) + J, M'First (2) + J);
|
|
Row := J;
|
|
for I in J .. Order - 1 loop
|
|
if M (M'First (1) + I, M'First (2) + J) > Max then
|
|
Max := M (M'First (1) + I, M'First (2) + J);
|
|
Row := I;
|
|
end if;
|
|
end loop;
|
|
if J /= Row then
|
|
-- swap rows J and Row
|
|
Swap_Rows (Result, J, Row);
|
|
end if;
|
|
end loop;
|
|
return Result;
|
|
end Pivoting_Matrix;
|
|
|
|
procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix) is
|
|
use type Matrix.Real_Matrix, Matrix.Real;
|
|
Order : constant Positive := A'Length (1);
|
|
A2 : Matrix.Real_Matrix (A'Range (1), A'Range (2));
|
|
S : Matrix.Real;
|
|
begin
|
|
L := (others => (others => 0.0));
|
|
U := (others => (others => 0.0));
|
|
P := Pivoting_Matrix (A);
|
|
A2 := P * A;
|
|
for J in 0 .. Order - 1 loop
|
|
L (L'First (1) + J, L'First (2) + J) := 1.0;
|
|
for I in 0 .. J loop
|
|
S := 0.0;
|
|
for K in 0 .. I - 1 loop
|
|
S := S + U (U'First (1) + K, U'First (2) + J) *
|
|
L (L'First (1) + I, L'First (2) + K);
|
|
end loop;
|
|
U (U'First (1) + I, U'First (2) + J) :=
|
|
A2 (A2'First (1) + I, A2'First (2) + J) - S;
|
|
end loop;
|
|
for I in J + 1 .. Order - 1 loop
|
|
S := 0.0;
|
|
for K in 0 .. J loop
|
|
S := S + U (U'First (1) + K, U'First (2) + J) *
|
|
L (L'First (1) + I, L'First (2) + K);
|
|
end loop;
|
|
L (L'First (1) + I, L'First (2) + J) :=
|
|
(A2 (A2'First (1) + I, A2'First (2) + J) - S) /
|
|
U (U'First (1) + J, U'First (2) + J);
|
|
end loop;
|
|
end loop;
|
|
end Decompose;
|
|
|
|
end Decomposition;
|