[Ada] Метод Якоби для решения системы линейных уравнений

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
 
with Ada.Numerics.Discrete_Random;
with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random;
 
with Ada.Numerics.Elementary_Functions; use  Ada.Numerics.Elementary_Functions;
 
procedure Main_Jacobi is
 
   Eps : constant Float := 0.00001;
 
   type Vector is array (Natural range <>) of Float;
   type Matrix is array (Natural range <>, Natural range <>) of Float;
 
   procedure Do_Jacobi (A : Matrix; F : Vector; X : in out Vector);
   function Is_Repeatable (A : Matrix) return Boolean;
   procedure Print (mx : Matrix);
   procedure Print (vec : Vector);
 
   procedure Do_Jacobi (A : Matrix; F : Vector; X : in out Vector) is
      TX : Vector (F'Range);
      Norm : Float;
   begin
      loop
         for i in TX'Range loop
            TX (i) := -F (i);
            for g in A'Range (2) loop
               if i /= g then
                  TX (i) := TX (i) + A (i, g) * X (g);
               end if;
            end loop;
            TX (i) := TX (i) / (-A (i, i));
         end loop;
 
         Norm := abs (X (0) - TX (0));
         for h in TX'Range loop
            if abs (X (h) - TX (h)) > Norm then
               Norm := abs (X (h) - TX (h));
            end if;
            X (h) := TX (h);
         end loop;
         exit when Norm <= Eps;
      end loop;
   end Do_Jacobi;
   
   function Is_Repeatable (A : Matrix) return Boolean is
      Ek : Float := 0.0;
   begin
      for i in A'Range (1) loop
         for j in A'Range (2) loop
            if i /= j then
               Ek := Ek + ((-A (i, j)) / A (i, i)) ** 2;
            end if;
         end loop;
      end loop;
      Ada.Float_Text_IO.Put (Ek, Fore => 8, Aft => 3, Exp => 0);
      Ada.Text_IO.Put_Line ("");
      Ek := Sqrt (Ek);
      Ada.Float_Text_IO.Put (Ek, Fore => 8, Aft => 3, Exp => 0);
      Ada.Text_IO.Put_Line ("");
      return (Ek < 1.0);
   end Is_Repeatable;
   
   procedure Print (mx : Matrix) is
   begin
      for i in mx'Range (1) loop
         for j in mx'Range (2) loop
            Ada.Float_Text_IO.Put (mx (i, j), Fore => 8, Aft => 3, Exp => 0);
         end loop;
         Ada.Text_IO.Put_Line ("");
      end loop;
   end Print;
 
   procedure Print (vec : Vector) is
   begin
      for i in vec'Range loop
         Ada.Float_Text_IO.Put (vec (i), Fore => 8, Aft => 3, Exp => 0);
      end loop;
      Ada.Text_IO.Put_Line ("");
   end Print;
   
   subtype Matrix_Size_Type is Integer range 5 .. 10;
   
   package Int_Rnd is new Ada.Numerics.Discrete_Random (Result_Subtype => Matrix_Size_Type);
   use Int_Rnd;
   I_Gen : Int_Rnd.Generator;
   F_Gen : Ada.Numerics.Float_Random.Generator;
 
begin
   Reset (I_Gen);  --  Integer generator
   Reset (F_Gen);  --  Float generator
   
   declare
      n : constant Integer := Int_Rnd.Random (I_Gen);
      --  4;
      subtype Index is Integer range 0 .. n - 1;
      
      mx : Matrix (Index, Index) := (others => (others => 0.0));
      --  ((15.0, 2.0, -1.0, -1.0),
      --  (1.0, -10.0, -1.0, -2.0),
      --  (2.0, 1.0, 12.0, 1.0),
      --  (1.0, 1.0, 1.0, 11.0));
 
      res : Vector (Index) := (others => 0.0);
      --  (22.0, -14.0, -10.0, -20.0);
      
      X : Vector (Index) := (others => 0.0);
 
   begin
      for i in Index'Range loop
         res (i) := 10.0 * (0.5 - Random (F_Gen));
         for j in Index'Range loop
            mx (i, j) := 10.0 * (0.5 - Random (F_Gen));
         end loop;
      end loop;
      
      Put_Line ("Матрица : ");
      Print (mx);
      Put_Line ("Вектор-столбец : ");
      Print (res);         
      
      if Is_Repeatable (mx) then   
         Do_Jacobi (mx, res, X);
         Put_Line ("Результат : ");
         Print (X);
      else
         Put_Line ("Ek > 1 : Итерационный процесс расходится");
      end if;
   end;
   
end Main_Jacobi;