10.4.2.9. Implementing the Generic Matrix Package

From the exploration of HB.Vectors and Show_Matrices, the package implementation is easy to understand.

     1 package body HB.Matrices_Generic is
     2
     3   function “*”
     4     (K : ElementType; Right : Vector) return Vector is
     5     Result: Vector(Right’Range);
     6   begin
     7     for I in Right’Range loop
     8       Result(I) := K * Right(I);
     9     end loop;
    10     return Result;
    11   end “*”;
    12
    13   function “+” (Left, Right : Vector) return Vector is
    14     Result : Vector(Left’Range);
    15   begin
    16     if Left’Length /= Right’Length then
    17       raise Bounds_Error;
    18     else
    19       for I in Left’Range loop
    20         Result(I) :=
    21           Left(I) + Right(I - Left’First + Right’First);
    22       end loop;
    23       return Result;
    24     end if;
    25   end “+”;
    26
    27   function “*”
    28     (Left, Right : Vector) return ElementType is
    29     Sum : ElementType := Zero;
    30   begin
    31     if Left’Length /= Right’Length then
    32       raise Bounds_Error;
    33     else
    34       for I in Left’Range loop
    35         Sum := Sum +
    36           Left(I) * Right(I - Left’First + Right’First);
    37       end loop;
    38       return Sum;
    39     end if;
    40   end “*”;
    41
    42   function “*”
    43     (K : ElementType; Right : Matrix) return Matrix is
    44     Result: Matrix(Right’Range(1), Right’Range(2));
    45   begin
    46     for I in Right’Range(1) loop
    47       for J in Right’Range(2) loop
    48         Result(I,J) := K * Right(I,J);
    49       end loop;
    50     end loop;
    51     return Result;
    52   end “*”;
    53
    54   function Transpose (Left: Matrix) return Matrix is
    55     Result: Matrix(Left’Range(2), Left’Range(1));
    56   begin
    57     for I in Result’Range(1) loop
    58       for J in Result’Range(2) loop
    59         Result(I,J) := Left(J,I);
    60       end loop;
    61     end loop;
    62     return Result;
    63   end Transpose;
    64
    65   function “+” (Left, Right : Matrix) return Matrix is
    66     Result : Matrix(Left’Range(1), Left’Range(2));
    67   begin
    68     if Left’Length(1) /= Right’Length(1) or
    69        Left’Length(2) /= Right’Length(2) then
    70       raise Bounds_Error;
    71     else
    72       for I in Left’Range(1) loop
    73         for J in Left’Range(2) loop
    74           Result(I, J) := Left(I, J) +
    75             Right(I - Left’First(1) + Right’First(1),
    76                   J - Left’First(2) + Right’First(2));
    77         end loop;
    78       end loop;
    79       return Result;
    80     end if;
    81   end “+”;
    82
    83   function “*” (Left, Right : Matrix) return Matrix is
    84     Result: Matrix(Left’Range(1), Right’Range(2)) :=
    85       (others => (others => Zero));
    86   begin
    87     if Left’Length(2) /= Right’Length(1) then
    88       raise Bounds_Error;
    89     else
    90       for RowL in Left’Range(1) loop
    91         for ColR in Right’Range(2) loop
    92           for ColL in Left’Range(2) loop
    93             Result(RowL, ColR) := Result(RowL, ColR) +
    94               Left(RowL, ColL) *
    95               Right((ColL-Left’First(2))+Right’First(1),
    96                     ColR);
    97           end loop;
    98         end loop;
    99       end loop;
   100     end if;
   101     return Result;
   102   end “*”;
   103
   104   function “*”
   105     (Left: Vector; Right: Matrix) return Vector is
   106     Result: Vector(Right’Range(2)) := (others => Zero);
   107   begin
   108     if Left’Length /= Right’Length(1) then
   109       raise Bounds_Error;
   110     else
   111       for C in Right’Range(2) loop
   112         for R in Right’Range(1) loop
   113           Result(C) := Result(C) + Right(R, C) *
   114             Left(Left’First + (R-Right’First(1)));
   115         end loop;
   116       end loop;
   117     end if;
   118     return Result;
   119   end “*”;
   120
   121   function “*”
   122     (Left: Matrix; Right: Vector) return Vector is
   123     Result: Vector(Left’Range(1)) := (others => Zero);
   124   begin
   125     if Left’Length(2) /= Right’Length then
   126       raise Bounds_Error;
   127     else
   128       for R in Left’Range(1) loop
   129         for C in Left’Range(2) loop
   130           Result(R) := Result(R) + Left(R, C) *
   131             Right(Right’First + (C-Left’First(2)));
   132         end loop;
   133       end loop;
   134     end if;
   135     return Result;
   136   end “*”;
   137
   138 end HB.Matrices_Generic;

I do not belabor the algorithms embodied in the various operations; these are well-known and straightforward algorithms from matrix algebra. From the Ada perspective, I point out the declaration (lines 84-85)

   Result: Matrix(Left’Range(1), Right’Range(2)) :=
             (others => (others => Zero));

which “sizes” a result matrix in terms of the Left and Right matrix operands and further initializes all elements of this matrix to 0, the parameter representing the zero of the element type. This aggregate initialization is very concise and can be implemented efficiently because the programmer leaves it up to the compiler just how all the elements are initialized at execution time.

Finally, lines 87-88 give the conformability condition for matrix multiplication: Left must have as many columns as Right has rows.