LU decomposition: Difference between revisions

=={{header|Ada}}==

decompose.ads:

<lang Ada>with Ada.Numerics.Generic_Real_Arrays;

generic

with package Matrix is new Ada.Numerics.Generic_Real_Arrays (<>);

package Decomposition is

-- decompose a square matrix A by PA = LU

procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix);

end Decomposition;</lang>

decompose.adb:

<lang 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;</lang>

Example usage:

<lang Ada>with Ada.Numerics.Real_Arrays;

with Ada.Text_IO;

with Decomposition;

procedure Decompose_Example is

package Real_Decomposition is new Decomposition

(Matrix => Ada.Numerics.Real_Arrays);

package Real_IO is new Ada.Text_IO.Float_IO (Float);

procedure Print (M : Ada.Numerics.Real_Arrays.Real_Matrix) is

begin

for Row in M'Range (1) loop

for Col in M'Range (2) loop

Real_IO.Put (M (Row, Col), 3, 2, 0);

end loop;

Ada.Text_IO.New_Line;

end loop;

end Print;

Example_1 : Ada.Numerics.Real_Arrays.Real_Matrix :=

((1.0, 3.0, 5.0),

(2.0, 4.0, 7.0),

(1.0, 1.0, 0.0));

P_1, L_1, U_1 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_1'Range (1),

Example_1'Range (2));

Example_2 : Ada.Numerics.Real_Arrays.Real_Matrix :=

((11.0, 9.0, 24.0, 2.0),

(1.0, 5.0, 2.0, 6.0),

(3.0, 17.0, 18.0, 1.0),

(2.0, 5.0, 7.0, 1.0));

P_2, L_2, U_2 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_2'Range (1),

Example_2'Range (2));

begin

Real_Decomposition.Decompose (A => Example_1,

P => P_1,

L => L_1,

U => U_1);

Real_Decomposition.Decompose (A => Example_2,

P => P_2,

L => L_2,

U => U_2);

Ada.Text_IO.Put_Line ("Example 1:");

Ada.Text_IO.Put_Line ("A:"); Print (Example_1);

Ada.Text_IO.Put_Line ("L:"); Print (L_1);

Ada.Text_IO.Put_Line ("U:"); Print (U_1);

Ada.Text_IO.Put_Line ("P:"); Print (P_1);

Ada.Text_IO.New_Line;

Ada.Text_IO.Put_Line ("Example 2:");

Ada.Text_IO.Put_Line ("A:"); Print (Example_2);

Ada.Text_IO.Put_Line ("L:"); Print (L_2);

Ada.Text_IO.Put_Line ("U:"); Print (U_2);

Ada.Text_IO.Put_Line ("P:"); Print (P_2);

end Decompose_Example;</lang>

Output:

<pre>Example 1:

A:

1.00 3.00 5.00

2.00 4.00 7.00

1.00 1.00 0.00

L:

1.00 0.00 0.00

0.50 1.00 0.00

0.50 -1.00 1.00

U:

2.00 4.00 7.00

0.00 1.00 1.50

0.00 0.00 -2.00

P:

0.00 1.00 0.00

1.00 0.00 0.00

0.00 0.00 1.00

Example 2:

A:

11.00 9.00 24.00 2.00

1.00 5.00 2.00 6.00

3.00 17.00 18.00 1.00

2.00 5.00 7.00 1.00

L:

1.00 0.00 0.00 0.00

0.27 1.00 0.00 0.00

0.09 0.29 1.00 0.00

0.18 0.23 0.00 1.00

U:

11.00 9.00 24.00 2.00

0.00 14.55 11.45 0.45

0.00 0.00 -3.47 5.69

0.00 0.00 0.00 0.51

P:

1.00 0.00 0.00 0.00

0.00 0.00 1.00 0.00

0.00 1.00 0.00 0.00

0.00 0.00 0.00 1.00</pre>