problem with Real_Matrix*Real_Matrix

problem with Real_Matrix*Real_Matrix

[reinkor]

Dear All,

I am confused about my result from the multiplication r*M where
r is Real_Vector and M is Real_Matrix. Enclosed is my test program
including comments.  Could someone look at it and maybe help me out
of my bubble?


I use gcc-ada-4.7-2.1.1.x86_64 under OpenSuse 12.2

reinert


with Text_IO;
use  Text_IO;
with   ada.numerics.Generic_Real_Arrays;

with Ada.Strings.Fixed;
use  Ada.Strings,Ada.Strings.Fixed;


procedure test3 is

  type Real is new Float;

  package Real_Io is new Text_IO.Float_Io   (Real);
  use Real_Io;

  package Int_Io is new Text_IO.Integer_Io (Integer);
  use Int_Io;

  package GRA is new Ada.Numerics.Generic_Real_Arrays(Real);
  use GRA;

  procedure Put (X : Real_Vector) is
  begin
      for I in X'Range (1) loop
          Put (X (I),6,1,0);
      end loop;
  end Put;

  M : Real_Matrix := ((1.0, 2.0),
                      (3.0, 4.0));

  r : Real_Vector := (10.0,1000.0);

begin

   Put("** r*M : ");New_Line;
   Put(r*M);

-- Extect: r*M = (10*1 + 1000*3, 10*2 + 1000*4) = (3010,4020)
-- i.e. should r*M consists of the inner product between r and 
-- the columns of M ?   But my program gives: 
--
-- ** r*M : 
--     40.0  6000.0
--      

end test3;

Re: problem with Real_Matrix*Real_Matrix

[Niklas Holsti]

On 12-09-17 16:37 , reinkor wrote:
> Dear All,
> 
> I am confused about my result from the multiplication r*M where
> r is Real_Vector and M is Real_Matrix. Enclosed is my test program
> including comments.  Could someone look at it and maybe help me out
> of my bubble?
> 
> 
> I use gcc-ada-4.7-2.1.1.x86_64 under OpenSuse 12.2

I get the expected result:

** r*M :
  3010.0  4020.0

on my Mac, with GNAT 4.5.0 20100221 (experimental) [trunk revision 156937].

I suspect your Ada installation is broken, but I don't have any advice
on how to fix it, sorry.

-- 
Niklas Holsti
Tidorum Ltd
niklas holsti tidorum fi
      .      @       .


> 
> reinert
> 
> 
> with Text_IO;
> use  Text_IO;
> with   ada.numerics.Generic_Real_Arrays;
> 
> with Ada.Strings.Fixed;
> use  Ada.Strings,Ada.Strings.Fixed;
> 
> 
> procedure test3 is
> 
>   type Real is new Float;
> 
>   package Real_Io is new Text_IO.Float_Io   (Real);
>   use Real_Io;
> 
>   package Int_Io is new Text_IO.Integer_Io (Integer);
>   use Int_Io;
> 
>   package GRA is new Ada.Numerics.Generic_Real_Arrays(Real);
>   use GRA;
> 
>   procedure Put (X : Real_Vector) is
>   begin
>       for I in X'Range (1) loop
>           Put (X (I),6,1,0);
>       end loop;
>   end Put;
> 
>   M : Real_Matrix := ((1.0, 2.0),
>                       (3.0, 4.0));
> 
>   r : Real_Vector := (10.0,1000.0);
> 
> begin
> 
>    Put("** r*M : ");New_Line;
>    Put(r*M);
> 
> -- Extect: r*M = (10*1 + 1000*3, 10*2 + 1000*4) = (3010,4020)
> -- i.e. should r*M consists of the inner product between r and 
> -- the columns of M ?   But my program gives: 
> --
> -- ** r*M : 
> --     40.0  6000.0
> --      
> 
> end test3;
> 
> 

Re: problem with Real_Matrix*Real_Matrix

[Anatoly Chernyshev]

Got the same mistake on Win7/GNAT 2012. As I'm doing matrix calculations too, I had to investigate, and here's the culprit from s-gearop.adb:
---------------------------------------------------------------
   function Vector_Matrix_Product
     (Left  : Left_Vector;
      Right : Matrix) return Result_Vector
   is
   begin
      return R : Result_Vector (Right'Range (2)) do
         if Left'Length /= Right'Length (2) then
            raise Constraint_Error with
              "incompatible dimensions in vector-matrix multiplication";
         end if;

         for J in Right'Range (2) loop
            declare
               S : Result_Scalar := Zero;

            begin
               for K in Right'Range (1) loop
                  S := S + Left (J - Right'First (1)
                                   + Left'First) * Right (K, J);
               end loop;

               R (J) := S;
            end;
         end loop;
      end return;
   end Vector_Matrix_Product;
---------------------------------------------------------------
S := S + Left (J - Right'First (1) + Left'First) * Right (K, J); - in this fragment Left is kept constant within the K-cycle. Replacing Left'First with K would correct the problem. 

I believe, it was straightened for some reason for the Mac distribution; would be interesting to compare.

Re: problem with Real_Matrix*Real_Matrix

[Simon Wright]

Anatoly Chernyshev <achernyshev@gmail.com> writes:

> Got the same mistake on Win7/GNAT 2012. As I'm doing matrix
> calculations too, I had to investigate, and here's the culprit from
> s-gearop.adb:
> ---------------------------------------------------------------
> function Vector_Matrix_Product (Left : Left_Vector; Right : Matrix)
> return Result_Vector is begin return R : Result_Vector (Right'Range
> (2)) do if Left'Length /= Right'Length (2) then raise Constraint_Error
> with "incompatible dimensions in vector-matrix multiplication"; end
> if;
>
>          for J in Right'Range (2) loop
>             declare
>                S : Result_Scalar := Zero;
>
>             begin
>                for K in Right'Range (1) loop
>                   S := S + Left (J - Right'First (1)
>                                    + Left'First) * Right (K, J);
>                end loop;
>
>                R (J) := S; end; end loop; end return; end
> Vector_Matrix_Product;
> --------------------------------------------------------------- 
> S := S + Left (J - Right'First (1) + Left'First) * Right (K, J); - in
> this fragment Left is kept constant within the K-cycle. Replacing
> Left'First with K would correct the problem.
>
> I believe, it was straightened for some reason for the Mac
> distribution; would be interesting to compare.

You're right about the location of the problem, but I think the fix is
instead to change

                  S := S + Left (J - Right'First (1)
to
                  S := S + Left (K - Right'First (1)

I've filed a bug in the GCC Bugzilla against GCC 4.7.0:
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=54614

Re: problem with Real_Matrix*Real_Matrix

[Anatoly Chernyshev]

On Tuesday, September 18, 2012 2:50:33 PM UTC+4, Simon Wright wrote:
> I think the fix is
> 
> instead to change
> 
> 
> 
>                   S := S + Left (J - Right'First (1)
> 
> to
> 
>                   S := S + Left (K - Right'First (1)

True.

Re: problem with Real_Matrix*Real_Matrix

[Egil Høvik]

On Monday, September 17, 2012 3:37:15 PM UTC+2, reinkor wrote:
> 
>    Put(r*M);
> 

As Anatoly pointed out, the bug is in Vector_Matrix_Product.
There's also a Matrix_Vector_Product, so you can probably work
around the bug by doing M*r instead.

-- 
~egilhh

Re: problem with Real_Matrix*Real_Matrix

[Anatoly Chernyshev]

On Tuesday, September 18, 2012 10:16:21 PM UTC+4, Egil Høvik wrote:
> On Monday, September 17, 2012 3:37:15 PM UTC+2, reinkor wrote:
> >    Put(r*M);
> 
> As Anatoly pointed out, the bug is in Vector_Matrix_Product.
> 
> There's also a Matrix_Vector_Product, so you can probably work
> 
> around the bug by doing M*r instead.

These operations are not commutative, yet if there's one bug, there must be another. For the time being I would suggest to switch to something tested by time, like lapack.

Re: problem with Real_Matrix*Real_Matrix

[Simon Wright]

Anatoly Chernyshev <achernyshev@gmail.com> writes:

> yet if there's one bug, there must be another.

Not sure that strictly follows! It does lead one to wonder about the
unit tests that either don't exist or were wrongly "fixed" to match the
changed behaviour.

> For the time being I would suggest to switch to something tested by
> time, like lapack.

In this case, that would be BLAS! (*gemv).

You can get this by using GCC 4.6 or 4.5, or GNAT GPL 2011.

Re: problem with Real_Matrix*Real_Matrix

[AdaMagica]

It's a rather scary bug, isn't it. In an operation like this!

Were AdaCore sleeping?

Re: problem with Real_Matrix*Real_Matrix

[reinkor]

On Tuesday, September 18, 2012 10:33:42 PM UTC+2, AdaMagica wrote:
> It's a rather scary bug, isn't it. In an operation like this!
> 
> 
> 
> Were AdaCore sleeping?

But they have some formalized QA? N (>>1) tests to pass through at
every release?

reinert

Re: problem with Real_Matrix*Real_Matrix

[Simon Wright]

reinkor <reinkor@gmail.com> writes:

> On Tuesday, September 18, 2012 10:33:42 PM UTC+2, AdaMagica wrote:
>> It's a rather scary bug, isn't it. In an operation like this!
>> 
>> Were AdaCore sleeping?
>
> But they have some formalized QA? N (>>1) tests to pass through at
> every release?

FSF GCC has the ACATS tests (not the latest version, though) + some
GNAT-specific tests (there are 947 source files, some tests require as
many as 3 sources).

AdaCore tell us that they have, in addition, test cases that they don't
make public (a lot of them are from bugs reported by customers and may
contain sensitive material); and that all these tests are run on a
regular basis.

I would have expected there to be a set of tests about numerics, and I'd
have expected these tests to be re-run after the refactoring that
happened at the GCC 4.6-4.7 change.

Re: problem with Real_Matrix*Real_Matrix

[Ken Thomas]

On Wednesday, 19 September 2012 09:17:17 UTC+1, Simon Wright  wrote:
> reinkor <reinkor@gmail.com> writes:
> 
> 
> 
> > On Tuesday, September 18, 2012 10:33:42 PM UTC+2, AdaMagica wrote:
> 
> >> It's a rather scary bug, isn't it. In an operation like this!
> 
> >> 
> 
> >> Were AdaCore sleeping?
> 
> >
> 
> > But they have some formalized QA? N (>>1) tests to pass through at
> 
> > every release?
> 
> 
> 
> FSF GCC has the ACATS tests (not the latest version, though) + some
> 
> GNAT-specific tests (there are 947 source files, some tests require as
> 
> many as 3 sources).
> 
> 
> 
> AdaCore tell us that they have, in addition, test cases that they don't
> 
> make public (a lot of them are from bugs reported by customers and may
> 
> contain sensitive material); and that all these tests are run on a
> 
> regular basis.
> 
> 
> 
> I would have expected there to be a set of tests about numerics, and I'd
> 
> have expected these tests to be re-run after the refactoring that
> 
> happened at the GCC 4.6-4.7 change.

This is indeed a pity. I have run some performance tests on a Windows machine (XP Intel 2 CPU 2GB ram). On multiplying two 500 by 500 matrices (Long_Float), a rate of 966 MFlops/Sec was achieved. This is good.