From mboxrd@z Thu Jan  1 00:00:00 1970
X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me
X-Spam-Level: 
X-Spam-Status: No, score=-0.3 required=5.0 tests=BAYES_00,
	REPLYTO_WITHOUT_TO_CC autolearn=no autolearn_force=no version=3.4.4
X-Google-Language: ENGLISH,ASCII-7-bit
X-Google-Thread: 103376,701270dccdaf7b2
X-Google-Attributes: gid103376,public
X-Google-ArrivalTime: 2003-06-01 02:16:31 PST
Path: 
 archiver1.google.com!news1.google.com!newsfeed.stanford.edu!news.uchicago.edu!newsfeed.cs.wisc.edu!144.212.100.101.MISMATCH!newsfeed!fu-berlin.de!uni-berlin.de!dialin-145-254-039-218.arcor-ip.NET!not-for-mail
From: "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de>
Newsgroups: comp.lang.ada
Subject: Re: In place matrix multiplication?
Date: Sun, 01 Jun 2003 11:18:59 +0200
Organization: At home
Message-ID: <bbcg9c$7u9v9$2@ID-77047.news.dfncis.de>
References: <1JcCa.16479$Mu3.361103@newsfep4-glfd.server.ntli.net>
Reply-To: mailbox@dmitry-kazakov.de
NNTP-Posting-Host: dialin-145-254-039-218.arcor-ip.net (145.254.39.218)
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7Bit
X-Trace: fu-berlin.de 1054458989 8333289 145.254.39.218 (16 [77047])
User-Agent: KNode/0.7.1
Xref: archiver1.google.com comp.lang.ada:38276
Date: 2003-06-01T11:18:59+02:00
List-Id: <comp.lang.ada>

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