Re: In place matrix multiplication?

["Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de>]

chris.danx wrote:

> Hi,
> 
> Does anyone know how to do this?  At first I thought it'd be something
> along the lines of...
> 
>     procedure Multiply (A : in out Matrix;
>                         B : in     Matrix) is
>        T : Scalar;
>        Z : Scalar := Create(0);
>     begin
>        for rCount in Matrix_Size'range loop
>           for cCount in Matrix_Size'range loop
>              T := Z;
>              for iCount in Matrix'range loop
>                 T := T + A.Mx(rCount, iCount) * B.Mx(iCount, cCount);
>              end loop;
>              A.Mx(rCount, iCount) := T;
>           end loop;
>        end loop;
>     end Multiply;
> 
> But this is obviously wrong.  Is there an algorithm to do this?  I
> googled for it, but haven't found anything yet.  I gather it's possible
> because some libraries do it.  They just don't say how they do so.

Something like this, I presume:

generic
   -- Some field of matrix elements:
   type Element is private;
   Zero : Element; -- Null of Element's field
   with function "+" (L, R : Element) return Element;
   with function "*" (L, R : Element) return Element;
   -- Some index type:
   type Index is (<>);
   -- Two-dimensional array over the field:
   type Matrix is array (Index range <>, Index range <>) of Element;
procedure Mul (A : in out Matrix; B : Matrix) is
begin
   if (  A'First (1) /= B'First (2)  -- A*B does not exist
      or else
         A'Last (1)  /= B'Last (2)
      or else
         A'First (1) /= A'First (2)  -- Can't store the result
      or else
         A'Last  (1) /= A'Last (2)
      )
   then
      raise Constraint_Error;
   end if;
   declare
      Temp : array (Index range A'Range (2));
      Item : Element;
   begin
      for Row in A'Range (1) loop
         -- Copy the row of A into Temp
         for Column in A'Range (2) loop
            Temp (Column) := A (Row, Column);
         end loop;
         -- Evaluate A*B (Row, *)
         for Column in A'Range (2) loop
            -- Evaluate A*B (Row, Column) in Item
            Item := Zero;
            for I in Temp'Range loop
               Item := Item + Temp (I) * B (I, Column);
            end loop;
            -- Store
            A (Row, Column) := Item;
         end loop;      
      end loop;
   end;
end Mul;

If you do not want to expose Zero, you can replace

-- Evaluate A*B (Row, Column) in Item
...

with a bit more cryptic:

Item := Temp (Temp'First) * B (Temp'First, Column);
if Temp'First /= Temp'Last then -- To avoid exception from 'Succ
   for I in Index'Succ (Temp'First)..Temp'Last loop
      Item := Item + Temp (I) * B (I, Column);
   end loop;
end if;

--
Regards,
Dmitry A. Kazakov
www.dmitry-kazakov.de