Source code in the ada programming language

generic
   type Element_Type is private;
   Zero : Element_Type;
   with function "-" (Left, Right : in Element_Type) return Element_Type is <>;
   with function "*" (Left, Right : in Element_Type) return Element_Type is <>;
   with function "/" (Left, Right : in Element_Type) return Element_Type is <>;
package Matrices is
   type Matrix is
     array (Positive range <>, Positive range <>) of Element_Type;
   function Reduced_Row_Echelon_form (Source : Matrix) return Matrix;
end Matrices;


package body Matrices is
   procedure Swap_Rows (From : in out Matrix; First, Second : in Positive) is
      Temporary : Element_Type;
   begin
      for Col in From'Range (2) loop
         Temporary          := From (First, Col);
         From (First, Col)  := From (Second, Col);
         From (Second, Col) := Temporary;
      end loop;
   end Swap_Rows;

   procedure Divide_Row
     (From    : in out Matrix;
      Row     : in Positive;
      Divisor : in Element_Type)
   is
   begin
      for Col in From'Range (2) loop
         From (Row, Col) := From (Row, Col) / Divisor;
      end loop;
   end Divide_Row;

   procedure Subtract_Rows
     (From                : in out Matrix;
      Subtrahend, Minuend : in Positive;
      Factor              : in Element_Type)
   is
   begin
      for Col in From'Range (2) loop
         From (Minuend, Col) := From (Minuend, Col) -
                                From (Subtrahend, Col) * Factor;
      end loop;
   end Subtract_Rows;

   function Reduced_Row_Echelon_form (Source : Matrix) return Matrix is
      Result : Matrix   := Source;
      Lead   : Positive := Result'First (2);
      I      : Positive;
   begin
      Rows : for Row in Result'Range (1) loop
         exit Rows when Lead > Result'Last (2);
         I := Row;
         while Result (I, Lead) = Zero loop
            I := I + 1;
            if I = Result'Last (1) then
               I    := Row;
               Lead := Lead + 1;
               exit Rows when Lead = Result'Last (2);
            end if;
         end loop;
         if I /= Row then
            Swap_Rows (From => Result, First => I, Second => Row);
         end if;
         Divide_Row
           (From    => Result,
            Row     => Row,
            Divisor => Result (Row, Lead));
         for Other_Row in Result'Range (1) loop
            if Other_Row /= Row then
               Subtract_Rows
                 (From       => Result,
                  Subtrahend => Row,
                  Minuend    => Other_Row,
                  Factor     => Result (Other_Row, Lead));
            end if;
         end loop;
         Lead := Lead + 1;
      end loop Rows;
      return Result;
   end Reduced_Row_Echelon_form;
end Matrices;


with Matrices;
with Ada.Text_IO;
procedure Main is
   package Float_IO is new Ada.Text_IO.Float_IO (Float);
   package Float_Matrices is new Matrices (
      Element_Type => Float,
      Zero => 0.0);
   procedure Print_Matrix (Matrix : in Float_Matrices.Matrix) is
   begin
      for Row in Matrix'Range (1) loop
         for Col in Matrix'Range (2) loop
            Float_IO.Put (Matrix (Row, Col), 0, 0, 0);
            Ada.Text_IO.Put (' ');
         end loop;
         Ada.Text_IO.New_Line;
      end loop;
   end Print_Matrix;
   My_Matrix : Float_Matrices.Matrix :=
     ((1.0, 2.0, -1.0, -4.0),
      (2.0, 3.0, -1.0, -11.0),
      (-2.0, 0.0, -3.0, 22.0));
   Reduced   : Float_Matrices.Matrix :=
      Float_Matrices.Reduced_Row_Echelon_form (My_Matrix);
begin
   Print_Matrix (My_Matrix);
   Ada.Text_IO.Put_Line ("reduced to:");
   Print_Matrix (Reduced);
end Main;