Re: Interfacing Ada with C

[Simon Wright <simon@pushface.org>]

[-- Attachment #1: Type: text/plain, Size: 666 bytes --]

Jeffrey Carter <spam.jrcarter.not@spam.not.acm.org> writes:

> On 08/05/2010 10:24 AM, Ada novice wrote:
>>
>> They are much more accurate and better results of course. I wonder why
>> Matlab tends to be give very "accurate" answers. Normally 15 digits of
>> precision is also used there.
>
> You requested 15 decimal digits; your results are zero to 15 decimal
> places. If you output in non-scientific notation you'd get zero;
> presumably that's what Matlab is doing.

The attached update uses fixed notation, looks rather better.

Note that the reason for the lack of precision in the Test16 results is
likely because the inputs are only specified to 6 digits.


[-- Attachment #2: New test program --]
[-- Type: text/plain, Size: 5011 bytes --]

--  This package is free software; you can redistribute it and/or
--  modify it under terms of the GNU General Public License as
--  published by the Free Software Foundation; either version 3, or
--  (at your option) any later version.  It is distributed in the
--  hope that it will be useful, but WITHOUT ANY WARRANTY; without
--  even the implied warranty of MERCHANTABILITY or FITNESS FOR A
--  PARTICULAR PURPOSE.
--
--  You should have received a copy of the GNU General Public License
--  along with this program; see the file COPYING3.  If not, see
--  <http://www.gnu.org/licenses/>.
--
--  Copyright Simon Wright <simon@pushface.org>

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Arrays.Extensions;

procedure Test_Extensions is

   subtype My_Float is Long_Float;
   package My_Float_IO is new Float_IO (My_Float);
   use My_Float_IO;

   package Real_Arrays
   is new Ada.Numerics.Generic_Real_Arrays (My_Float);
   package Complex_Types
   is new Ada.Numerics.Generic_Complex_Types (My_Float);
   package Complex_Arrays
   is new Ada.Numerics.Generic_Complex_Arrays (Real_Arrays, Complex_Types);
   package Extensions
   is new Complex_Arrays.Extensions;

   use Real_Arrays;
   use Complex_Types;
   use Complex_Arrays;

begin

   declare
   --  Values in yc's example
      Input : constant Complex_Matrix
        := (((8.0, 0.0), (-1.0, 0.0), (-5.0, 0.0)),
            ((-4.0, 0.0), (4.0, 0.0), (-2.0, 0.0)),
            ((18.0, 0.0), (-5.0, 0.0), (-7.0, 0.0)));
      Result : Complex_Vector (1 .. Input'Length (1));
   begin

      Put_Line ("Values from <143ef70b-7e74-426b-a621-a5fd157849be"
                  & "@x21g2000yqa.googlegroups.com>");

      Result := Extensions.Eigenvalues (Input);

      for J in Result'Range loop
         Put (Result (J).Re, Exp => 0);
         Put (" ");
         Put (Result (J).Im, Exp => 0);
         New_Line;
      end loop;

      New_Line;

   end;

   declare
   --  Values in Test 16 of
   --  http://people.sc.fsu.edu/~jburkardt/f_src/lapack/lapack_OSX_prb_output.txt
      Z : constant Complex := (0.0, 0.0);
      A : constant Complex := (2.44949, 0.0);
      B : constant Complex := (3.16228, 0.0);
      C : constant Complex := (3.46410, 0.0);
      Input : constant Complex_Matrix
        := (
            1 => (Z, A, Z, Z, Z, Z, Z),
            2 => (A, Z, B, Z, Z, Z, Z),
            3 => (Z, B, Z, C, Z, Z, Z),
            4 => (Z, Z, C, Z, C, Z, Z),
            5 => (Z, Z, Z, C, Z, B, Z),
            6 => (Z, Z, Z, Z, B, Z, A),
            7 => (Z, Z, Z, Z, Z, A, Z)
           );
   begin

      Put_Line ("Values in Test16 of "
                  & "http://people.sc.fsu.edu/~jburkardt/f_src/lapack"
                  & "/lapack_OSX_prb_output.txt");

      -- GNAT: Eigenvalues of symmetrix complex matrix are real
      Put_Line ("using Complex_Arrays.Eigenvalues");
      declare
         Result : constant Real_Vector := Complex_Arrays.Eigenvalues (Input);
      begin
         for J in Result'Range loop
            Put (Result (J), Exp => 0);
            New_Line;
         end loop;
      end;
      New_Line;

      --  Extension: Eigenvalues of general complex matrix are complex.
      Put_Line ("using Extensions.Eigenvalues");
      declare
         Result : constant Complex_Vector := Extensions.Eigenvalues (Input);
      begin
         for J in Result'Range loop
            Put (Result (J).Re, Exp => 0);
            Put (" ");
            Put (Result (J).Im, Exp => 0);
            New_Line;
         end loop;
      end;
      New_Line;

   end;

   declare
   --  Values from http://en.wikipedia.org/wiki/Skew-symmetric_matrix
      Input : constant Complex_Matrix
        := (((0.0, 0.0), (2.0, 0.0), (-1.0, 0.0)),
            ((-2.0, 0.0), (0.0, 0.0), (-4.0, 0.0)),
            ((1.0, 0.0), (4.0, 0.0), (0.0, 0.0)));
      Result : Complex_Vector (1 .. Input'Length (1));
   begin

      Put_Line
        ("Values from http://en.wikipedia.org/wiki/Skew-symmetric_matrix");

      Result := Extensions.Eigenvalues (Input);

      for J in Result'Range loop
         Put (Result (J).Re, Exp => 0);
         Put (" ");
         Put (Result (J).Im, Exp => 0);
         New_Line;
      end loop;
      New_Line;

   end;

   declare
   --  Values from http://en.wikipedia.org/wiki/Orthogonal_matrix
      Input : constant Complex_Matrix
        := (((0.0, 0.0), (-0.8, 0.0), (-0.6, 0.0)),
            ((0.8, 0.0), (-0.36, 0.0), (0.48, 0.0)),
            ((0.6, 0.0), (0.48, 0.0), (-0.64, 0.0)));
      Result : Complex_Vector (1 .. Input'Length (1));
   begin

      Put_Line
        ("Values from http://en.wikipedia.org/wiki/Orthogonal_matrix");

      Result := Extensions.Eigenvalues (Input);

      for J in Result'Range loop
         Put (Result (J).Re, Exp => 0);
         Put (" ");
         Put (Result (J).Im, Exp => 0);
         New_Line;
      end loop;
      New_Line;

   end;

end Test_Extensions;