Метод прогонки — на языке Ада ]:->

-- WARNING: Run only on MS Windows (version 5.0 or older)!
--          Запускать только под MS Windows (версии 5.0 или старше)!
 
with Ada.Text_IO;                       use Ada.Text_IO;
with Ada.Numerics.Float_Random;
with Ada.Integer_Text_IO;
with Ada.Float_Text_IO;
with Interfaces.C;
with Ada.Numerics.Elementary_Functions;
 
--Программа находит решение x системы уранений Ax=f, где
--A - трёхдиагональная матрица, а f - вектор-столбец.
--Общий вид расширенной матрицы системы:
--c1 b1 0  ... 0 |f1
--a2 c2 b2 ... 0 |f2
--0  a3 c3 ... 0 |f3
--.  .  .  ... . |.
--.  .  .  ... . |.
--0  0  ... an cn|fn
--Поиск решения осуществляется методом прогонки.
 
procedure Main is
   type Vector is array (Integer range <>) of Float;
   Size       : Integer;
   Min_Random : Float := -2.0;
   Max_Random : Float := 2.0;
 
   procedure System_Command (s : in String) is
      -- выолнение командной строки
      use Interfaces.C;
      function Sys (Arg : char_array) return Integer;
      pragma Import (C, Sys, "system");
      Ret_Val : Integer;
   begin
      Ret_Val := Sys (To_C (s));
   end System_Command;
 
   procedure Set_Console_Encoding (v : in Interfaces.C.unsigned) is
      -- изменение кодировки консоли
      use Interfaces.C;
      function SetCons
        (Arg  : Interfaces.C.unsigned)
         return Interfaces.C.unsigned_short;
      pragma Import (Stdcall, SetCons, "SetConsoleCP");
      function SetConsOut
        (Arg  : Interfaces.C.unsigned)
         return Interfaces.C.unsigned_short;
      pragma Import (Stdcall, SetConsOut, "SetConsoleOutputCP");
      Ret_Val : unsigned_short;
   begin
      Ret_Val := SetCons (v);
      Ret_Val := SetConsOut (v);
   end Set_Console_Encoding;
 
   function First_Digit (v : Float) return Integer is
      -- поиск первой значащей цифры
      use Ada.Numerics.Elementary_Functions;
   begin
      if v = 0.0 then
         return Integer'First;
      end if;
      return Integer (Float'Truncation (Log (X => v, Base => 10.0)));
   end First_Digit;
 
   procedure Print_Vector
     (v              : in Vector;
      Print_New_Line : in Boolean := False)
   is
      -- вывод вектора на экран
      fzs  : Integer :=
        First_Digit (Float'Max (abs Max_Random * 4.0, abs Min_Random * 4.0));
      fore : Integer := (if Fzs>0 then Fzs+2 else 2);
      aft  : Integer := (if Fzs<0 then 3-Fzs else 3);
   begin
      Put ("(");
      for I in v'First .. v'Last - 1 loop
         Ada.Float_Text_IO.Put (v (I), fore, aft, 0);
         Put (",");
      end loop;
      Ada.Float_Text_IO.Put (v (v'Last), fore, aft, 0);
      if Print_New_Line then
         Put_Line (")");
      else
         Put (")");
      end if;
   end Print_Vector;
 
   procedure Print_Vector
     (v              : in Vector;
      Fzs            : in Integer;
      Print_New_Line : in Boolean := False;
      Inc_Aft        : in Boolean := False)
   is
      -- вывод вектора на экран необычным образом
      fore : Integer := (if Fzs>0 then Fzs+2 else 2);
      aft  : Integer := (if Fzs<0 then 1-Fzs else 0);
   begin
      if Inc_Aft then
         aft := aft + 3 - (if Fzs<0 then 1 else 0);
      end if;
      Put ("(");
      for I in v'First .. v'Last - 1 loop
         Ada.Float_Text_IO.Put
           (v (I),
            fore + (if v(I)<0.0 then 1 else 0),
            aft,
            0);
         Put (",");
      end loop;
      Ada.Float_Text_IO.Put
        (v (v'Last),
         fore + (if v(v'Last)<0.0 then 1 else 0),
         aft,
         0);
      if Print_New_Line then
         Put_Line (")");
      else
         Put (")");
      end if;
   end Print_Vector;
 
   procedure Random_Vector
     (v   : in out Vector;
      min : Float := -2.0;
      max : Float := 2.0)
   is
      -- создание случайного вектора
      use Ada.Numerics.Float_Random;
      seed : Generator;
   begin
      Reset (seed);
      delay 0.1;
      for I in v'Range loop
         v (I) := Random (seed) * (max - min) + min;
      end loop;
   end Random_Vector;
 
   function Read_Vector
     (help  : in String;
      n     : in Integer;
      first : in Integer := 1)
      return  Vector
   is
      -- чтение вектора из ввода пользователя
      v : Vector (first .. n + first - 1);
   begin
      Put_Line ("Введите вектор " & help & ":");
      for i in v'Range loop
         Put ("  " & help & "[" & i'Img & "]=");
         Ada.Float_Text_IO.Get (v (i));
         if v(i)/4.0>Max_Random then
            Max_Random:=v(i)/4.0;
         end if;
         if v(i)/4.0<Min_Random then
            Min_Random:=v(i)/4.0;
         end if;
      end loop;
      return v;
   end Read_Vector;
 
   procedure Print_3diagonal (a, b, c, f : in Vector) is
      -- вывод расширенной матрицы системы в консоль
      tmp, nula : Vector (c'Range) := (others => 0.0);
      asd       : Integer          := a'First - c'First;
      bsd       : Integer          := b'First - c'First;
      fsd       : Integer          := f'First - c'First;
   begin
      tmp (c'First .. c'First + 1) := Vector'(c (c'First), b (b'First));
      Print_Vector (tmp);
      Put ("|");
      Ada.Float_Text_IO.Put (f (f'First), 2, 3, 0);
      New_Line;
      for i in c'First + 1 .. c'Last - 1 loop
         tmp                  := nula;
         tmp (i - 1 .. i + 1) := Vector'(a (i + asd), c (i), b (i + bsd));
         Print_Vector (tmp);
         Put ("|");
         Ada.Float_Text_IO.Put (f (i + fsd), 2, 3, 0);
         New_Line;
      end loop;
      tmp                        := nula;
      tmp (c'Last - 1 .. c'Last) := (a (a'Last), c (c'Last));
      Print_Vector (tmp);
      Put ("|");
      Ada.Float_Text_IO.Put (f (f'Last), 2, 3, 0);
      New_Line;
   end Print_3diagonal;
 
   function Diagonal_Dominant (a, b, c : in Vector) return Boolean is
      -- проверка на диагональное преобладание
      asd : Integer := a'First - c'First;
      bsd : Integer := b'First - c'First;
   begin
      if (abs c (c'First) < abs b (b'First)) or
         (abs c (c'Last) < abs a (a'Last))
      then
         return False;
      end if;
      for i in c'First + 1 .. c'Last - 1 loop
         if abs c (i) < (abs b (i + bsd) + abs a (i + asd)) then
            return False;
         end if;
      end loop;
      return True;
   end Diagonal_Dominant;
 
   function Solve_3diagonal (a, b, c, f : in Vector) return Vector is
      -- собственно решение
      m   : Float;
      x   : Vector (f'Range);
      tc  : Vector (c'Range) := c;
      tf  : Vector (f'Range) := f;
      asd : Integer          := a'First - x'First - 1;
      bsd : Integer          := b'First - x'First;
      csd : Integer          := c'First - x'First;
   begin
      for i in x'First + 1 .. x'Last loop
         m      := a (i + asd) / tc (i - 1 + csd);
         tc (i) := tc (i + csd) - m * b (i - 1 + bsd);
         tf (i) := tf (i) - m * tf (i - 1);
      end loop;
      x (x'Last) := tf (tf'Last) / tc (tc'Last);
      for i in reverse x'First .. x'Last - 1 loop
         x (i) := (tf (i) - b (i + bsd) * x (i + 1)) / tc (i + csd);
      end loop;
      return x;
   end Solve_3diagonal;
 
   function Residual (a, b, c, f, x : in Vector) return Vector is
      -- вычисление невязки
      v   : Vector (f'Range) := (others => 0.0);
      asd : Integer          := a'First - f'First - 1;
      bsd : Integer          := b'First - f'First;
      csd : Integer          := c'First - f'First;
      xsd : Integer          := x'First - f'First;
   begin
      v (v'First) := c (c'First) * x (x'First) +
                     b (b'First) * x (x'First + 1) -
                     f (f'First);
      for i in f'First + 1 .. f'Last - 1 loop
         v (i) := a (i + asd) * x (i + xsd - 1) +
                  c (i + csd) * x (i + xsd) +
                  b (i + bsd) * x (i + xsd + 1) -
                  f (i);
      end loop;
      v (v'Last) := c (c'Last) * x (x'Last) +
                    a (a'Last) * x (x'Last - 1) -
                    f (f'Last);
      return v;
   end Residual;
 
   function Max_Value_Sim (v : in Vector) return Integer is
      -- поиск первой значащей цифры максимального элемента вектора
      tmp : Float;
   begin
      tmp := abs (v (v'First));
      for I in v'First + 1 .. v'Last loop
         tmp := Float'Max (abs (v (I)), tmp);
      end loop;
      return First_Digit (tmp);
   end Max_Value_Sim;
 
begin
   Set_Console_Encoding (1251);-- укажите здесь кодировку файла'
   Put_Line ("                         Метод прогонки на языке Ада");
   Put ("Введите число уравнений: ");
   Ada.Integer_Text_IO.Get (Size);
   if Size <= 1 then -- система из одного уравнения - это уж слишком
      Size := 2;
   end if;
   declare
      a : Vector (2 .. Size);
      -- вектор, описывающий диагональ, лежащую под главной
      b : Vector (1 .. Size - 1);
      -- вектор, описывающий диагональ, лежащую над главной
      c, f, x : Vector (1 .. Size);
      -- главная диагональ, правые части и решение
      tmp     : Character;
      Rnd_Nmb : Boolean;
   begin
      Put ("Задать значения вручную (д/н)? ");
      Get_Immediate (tmp);
      Rnd_Nmb :=
        (case tmp is
        when 'y'|'Y'|'д'|'Д'=>false,
        when others => true);
      if Rnd_Nmb then
         Put_Line ("Нет");
         Put ("min/4=");
         Ada.Float_Text_IO.Get (Min_Random);
         Put ("max/4=");
         Ada.Float_Text_IO.Get (Max_Random);
         --В главной диагонали и правой части значения из R(min,max),
         --остальные числа из R(min/4,max/4). Это сделано для увеличения
         --вероятности наличия строгого диагонального преобладания
         loop
            Random_Vector (a, Min_Random, Max_Random);
            Random_Vector (b, Min_Random, Max_Random);
            Random_Vector (c, Min_Random * 4.0, Max_Random * 4.0);
            exit when Diagonal_Dominant (a, b, c);
         end loop;
         Random_Vector (f, Min_Random * 4.0, Max_Random * 4.0);
      else
         Put_Line ("Да");
         loop
            a := Read_Vector ("a", a'Length, a'First);
            b := Read_Vector ("b", b'Length);
            c := Read_Vector ("c", c'Length);
            exit when Diagonal_Dominant (a, b, c);
            Put_Line ("Требуется строгое диагональное преобладание!");
            Put_Line ("Попробуйте ещё раз");
         end loop;
         f := Read_Vector ("f", f'Length);
      end if;
      Print_3diagonal (a, b, c, f);
      New_Line;
      x := Solve_3diagonal (a, b, c, f);
      Put ("Решение = ");
      Print_Vector (x, Max_Value_Sim (x), True, True);
      x := Residual (a, b, c, f, x);
      Put ("Невязка = ");
      Print_Vector (x, Max_Value_Sim (x), True);
      --если здесь сильно не нули, то либо система неустойчива, либо нет
      --диагонального преобладания, либо числа совсем уж плохие
   end;
   New_Line;
   Put_Line ("Выполнение программы завершено");
   System_Command ("pause");
end Main;