This is something of a newbie question..

I'm working with a two dimensional array of floating point values. Lets call
it A. The array has type

type Matrix is array(Positive range <>, Positive range <>) of Floating_Type;

I need to exchange two rows in this array. What I'd like to do is something
along these lines:

   Temp_Array      := A(I, 1 .. Size);
   A(I, 1 .. Size) := A(K, 1 .. Size);
   A(K, 1 .. Size) := Temp_Array;

The compiler (GNAT GPL 2009) has a problem with this syntax and, after looking
into it some, I think that's because slicing only works for one dimensional
arrays. Fair enough.

So I thought, "Perhaps A needs to be an array of arrays."

type Matrix is array(Positive range <>) of WHAT_EXACTLY?

Apparently the component type of an array needs to be fully constrained (which
again makes sense) yet I don't know the size I'll want to use at the point
where this type is declared.

So now I'm thinking that I'll have to write a procedure to explicity swap each
row element one at a time. Of course this is not a terrible thing, but I'm
wondering if there is a more elegant way that I'm missing. I have some
confidence that the compiler can optimize slice operations reasonably well.
I'm less confident about its ability to optimize element by element
operations (maybe I'm overly pessimistic).

Peter

Peter C. Chapin wrote:
> 
> type Matrix is array(Positive range <>, Positive range <>) of Floating_Type;
> 
> I need to exchange two rows in this array. What I'd like to do is something
> along these lines:
> 
>    Temp_Array      := A(I, 1 .. Size);
>    A(I, 1 .. Size) := A(K, 1 .. Size);
>    A(K, 1 .. Size) := Temp_Array;
> 
> The compiler (GNAT GPL 2009) has a problem with this syntax and, after looking
> into it some, I think that's because slicing only works for one dimensional
> arrays. Fair enough.

That's correct: slicing is only defined for one-D arrays.

> So I thought, "Perhaps A needs to be an array of arrays."
> 
> type Matrix is array(Positive range <>) of WHAT_EXACTLY?

Row_Vector?

> Apparently the component type of an array needs to be fully constrained (which
> again makes sense) yet I don't know the size I'll want to use at the point
> where this type is declared.

Correct again: array components must be definite.

I'm not aware of a solution other than the brute force approach.

-- 
Jeff Carter
"Gentlemen, you can't fight in here. This is the War Room!"
Dr. Strangelove
30

Peter C. Chapin wrote:
> This is something of a newbie question..
> 
> I'm working with a two dimensional array of floating point values. Lets call
> it A. The array has type
> 
> type Matrix is array(Positive range <>, Positive range <>) of Floating_Type;
> 
> I need to exchange two rows in this array. What I'd like to do is something
> along these lines:
> 
>    Temp_Array      := A(I, 1 .. Size);
>    A(I, 1 .. Size) := A(K, 1 .. Size);
>    A(K, 1 .. Size) := Temp_Array;
> 
> The compiler (GNAT GPL 2009) has a problem with this syntax and, after looking
> into it some, I think that's because slicing only works for one dimensional
> arrays. Fair enough.
> 
> So I thought, "Perhaps A needs to be an array of arrays."
> 
> type Matrix is array(Positive range <>) of WHAT_EXACTLY?
> 
> Apparently the component type of an array needs to be fully constrained (which
> again makes sense) yet I don't know the size I'll want to use at the point
> where this type is declared.

One solution -- which may not match your design in other respects -- is 
to make the Matrix an array of accesses to rows:

    type Row is array (Positive range <>) of Float;
    type Row_Ref is access Row;
    type Matrix is array (Positive range <>) of Row_Ref;

Exchanging two rows is then very quick:

    Temp_Ref := A(I);
    A(I)     := A(K);
    A(K)     := Temp_Ref;

However, you now have to allocate the rows using "new Row (1 .. Size)" 
and perhaps later deallocate them using Unchecked_Deallocation.

Another solution -- which perhaps you have considered -- is to use a 
generic package to define row and matrix types:

    generic Size : Positive;
    package Matrices
    is
       type Row is array (1 .. Size) of Float;
       type Matrix is array (1 .. Size) of Row;
    end Matrices;

However, this will probably force some other parts of your code to 
become generic, too, parametrized either by a Size, or by an instance of 
Matrices.

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

On Sun, 31 Jan 2010 21:11:12 -0500, Peter C. Chapin wrote:

> So I thought, "Perhaps A needs to be an array of arrays."
> 
> type Matrix is array(Positive range <>) of WHAT_EXACTLY?
> 
> Apparently the component type of an array needs to be fully constrained (which
> again makes sense)

Not really. Arrays should have been allowed to have discriminants. That was
missed in the language design. If they had discriminants you could write:

type Row is array (Positive range <>) of Float;
type Matrix (Width : Positive) is
   array (Positive range <>) of Row (1..Width);

(Of course  the array bounds should have been discriminants rather than
attributes)

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

I've never understood why Ada does not allow slicing in
multidimensional arrays. What are the safety issues involved? And how
is it safe to force the programmer into ad hoc methods?

Jerry

Niklas Holsti wrote:

> Another solution -- which perhaps you have considered -- is to use a 
> generic package to define row and matrix types:
> 
>     generic Size : Positive;
>     package Matrices
>     is
>        type Row is array (1 .. Size) of Float;
>        type Matrix is array (1 .. Size) of Row;
>     end Matrices;
> 
> However, this will probably force some other parts of your code to 
> become generic, too, parametrized either by a Size, or by an instance of 
> Matrices.
> 

Actually I hadn't considered this, but such an approach might work well for
me. I'll have to consider it more.

Thanks!

Peter

I presume this was because the elements in question are not necessarily 
adjacent. As such the code ends up being exactly the same (or sometimes 
worse) than an explicitly coded loop.

And honestly, I don't see any real advantage to slices in terms of safety 
(given that indexes are checked for in-range, as they are in Ada). I make 
just as many mistakes with slices as I do with loops -- I find it hard to 
figure out the ends of the two slices for assignment (which have to be the 
same length - mine only are about 1/2 of the time); it's about the same 
difficulty as figuring out the offset amount when writing a loop to copy 
items. So I see it as a wash, other than that a one-dimensional slice can 
generate much better code (using a direct memory-memory move the for entire 
amount). The safety win comes for whole object assignments, not parts.

                           Randy.

"Jerry" <lanceboyle@qwest.net> wrote in message 
news:ed36036c-8318-4f27-aaae-5329a8bfc83d@t31g2000prh.googlegroups.com...
> I've never understood why Ada does not allow slicing in
> multidimensional arrays. What are the safety issues involved? And how
> is it safe to force the programmer into ad hoc methods?
>
> Jerry 

"Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> wrote in message 
news:178rg3rch8qdu$.13cgkywb09x3p.dlg@40tude.net...
> On Sun, 31 Jan 2010 21:11:12 -0500, Peter C. Chapin wrote:
>
>> So I thought, "Perhaps A needs to be an array of arrays."
>>
>> type Matrix is array(Positive range <>) of WHAT_EXACTLY?
>>
>> Apparently the component type of an array needs to be fully constrained 
>> (which
>> again makes sense)
>
> Not really. Arrays should have been allowed to have discriminants. That 
> was
> missed in the language design.

Early Ada 9x had discriminants for arrays. The capability got dropped when 
lots of "nice-to-have" features got dropped from Ada 95. Such "array" types 
have to be implemented as discriminant dependent records anyway; there is no 
real hope of performance improvement from them, but a lot of implementation 
complication. (I recall I was one of the stronger opponents of this feature, 
mostly for implementation cost/benefit reasons.)

So it is completely wrong to say that "this was missed in the language 
design". The correct statement is that "this was explicitly rejected in the 
language design".

                                 Randy.

Jerry a �crit :
> I've never understood why Ada does not allow slicing in
> multidimensional arrays. What are the safety issues involved? And how
> is it safe to force the programmer into ad hoc methods?
> 
One-dimensional slices are simple and efficient. Multidimensional slices
are a can of worms.

I guess you are thinking about rectangular slices. But why stop there? A
slice comprising the main diagonal and the diagonal above and below can
be very useful for some calculations. Or a slice which is the triangular
part of the upper half...

AFAICT, Fortran-99 does provide this - and the syntax is so complicated
that nobody uses it. And implementation is also a nightmare.

When designing a programming language, you have to stop at some point.
The ratio (cost of implementation) / usefulness is a good measure for
this. I think the ratio was simply to high for this feature.

-- 
---------------------------------------------------------
           J-P. Rosen (rosen@adalog.fr)
Visit Adalog's web site at http://www.adalog.fr

On Feb 2, 1:52=A0am, Jean-Pierre Rosen <ro...@adalog.fr> wrote:
> Jerry a =E9crit :> I've never understood why Ada does not allow slicing i=
n
> > multidimensional arrays. What are the safety issues involved? And how
> > is it safe to force the programmer into ad hoc methods?
>
> One-dimensional slices are simple and efficient. Multidimensional slices
> are a can of worms.
>
> I guess you are thinking about rectangular slices. But why stop there? A
> slice comprising the main diagonal and the diagonal above and below can
> be very useful for some calculations. Or a slice which is the triangular
> part of the upper half...
>
> AFAICT, Fortran-99 does provide this - and the syntax is so complicated
> that nobody uses it. And implementation is also a nightmare.
>
> When designing a programming language, you have to stop at some point.
> The ratio (cost of implementation) / usefulness is a good measure for
> this. I think the ratio was simply to high for this feature.
>
> --
> ---------------------------------------------------------
> =A0 =A0 =A0 =A0 =A0 =A0J-P. Rosen (ro...@adalog.fr)
> Visit Adalog's web site athttp://www.adalog.fr

Well, yes, I was thinking of rectangular slices. No doubt the (cost of
implementation) / usefulness is high (and usage difficult) for non-
rectangular parts, but that is far less common than rectangular parts.
Python, Matlab/Octave, Igor Pro... all pull it off without too much
hassle (although Python asks you to imagine addressing the array by
the "cracks" between elements, as I recall--probably a disease of C-
style counting).

Jerry

On Feb 2, 2:23=A0pm, Jerry <lancebo...@qwest.net> wrote:
> On Feb 2, 1:52=A0am, Jean-Pierre Rosen <ro...@adalog.fr> wrote:
>
>
>
>
>
> > Jerry a =E9crit :> I've never understood why Ada does not allow slicing=
 in
> > > multidimensional arrays. What are the safety issues involved? And how
> > > is it safe to force the programmer into ad hoc methods?
>
> > One-dimensional slices are simple and efficient. Multidimensional slice=
s
> > are a can of worms.

> Well, yes, I was thinking of rectangular slices. No doubt the (cost of
> implementation) / usefulness is high (and usage difficult) for non-
> rectangular parts, but that is far less common than rectangular parts.
> Python, Matlab/Octave, Igor Pro... all pull it off without too much
> hassle (although Python asks you to imagine addressing the array by
> the "cracks" between elements, as I recall--probably a disease of C-
> style counting).

I don't know anything about any of those languages.  What I know about
Ada is that when multi-dimensional arrays are passed as parameters to
procedures, they're normally passed by reference (if they're large
enough).  Some arrays are required to be passed by reference; in other
cases, it's unspecified, but I'd expect any array to be passed by
reference unless it's a packed array of 32 Booleans or something.  So
say you have a procedure:

   type Matrix is array (Natural range <>, Natural range <>) of
Float;
   procedure Operation_On_Matrix (M : in out Matrix) is ...

You later declare a matrix:

   X : Matrix (1..10, 1..10);

and want to call the operation on a rectangular slice.  If the "slice"
consists of an entire row, or one or more consecutive rows:

   Operation_On_Matrix (X (3..4, 1..10));

this could be done as efficiently as if all of X were passed, since
the procedure would see it as an array of 20 consecutive floats (in a
typical implementation).  However, a slice consisting of one or more
columns:

   Operation_On_Matrix (X (1..10, 5..6));

or some smaller rectangle:

   Operation_On_Matrix (X (4..6, 4..6));

would be tricky, since the procedure now has to be told that there are
gaps between the "rows" of the array that it's seeing as a parameter.
This means making Operation_On_Matrix less efficient, since it has to
be given more information about each Matrix that comes in; and the
inefficiency has to be put there even if no caller ever used a 2-D
rectangular slice, ever.  Maybe the efficiency hit isn't all that
large---I don't know.

It wouldn't really work to fix the language to say that rectangular
slices are allowed only if they include the entire index range in all
dimensions but the first, since that would fail for array types
declared with the Fortran convention.

I realize that the original question was about using assignment to
copy slices.  I suppose that, theoretically, the language could be
changed to allow multi-dimensional slices in assignments but not as
subprogram parameters.  Yuck.  I wouldn't want to add that kind of
inconsistency to the language.

Anyway, I don't know how serious these issues are, but this seems to
me to be a possible reason why adding this feature isn't as simple as
it sounds.

                                  -- Adam

On 1 f=E9v, 03:11, "Peter C. Chapin" <pcc482...@gmail.com> wrote:
> Apparently the component type of an array needs to be fully constrained (=
which
> again makes sense) yet I don't know the size
That's what generics are made for.
But you will have a single matrix type of a fixed size.

I you don't want a fixed size type, the best I will think about, would
be to create an abstract matrix type which will be implemented on a
one dimensional array. You could use slice in the implementation.

Ex. for a 2x2 matrix, you will use a 1 .. 4 array of Float. You may
define a swap row method which will move items from 1 .. 2 to 3 .. 4
and vice-versa.

You will be able to use slice optimized operations.

Tell me if you need a more concrete example (I confess I'm a bit terse
here).

On 1 f=E9v, 07:55, Niklas Holsti <niklas.hol...@tidorum.invalid> wrote:
> Another solution -- which perhaps you have considered -- is to use a
> generic package to define row and matrix types:
>
> =A0 =A0 generic Size : Positive;
> =A0 =A0 package Matrices
> =A0 =A0 is
> =A0 =A0 =A0 =A0type Row is array (1 .. Size) of Float;
> =A0 =A0 =A0 =A0type Matrix is array (1 .. Size) of Row;
> =A0 =A0 end Matrices;
>
> However, this will probably force some other parts of your code to
> become generic, too, parametrized either by a Size, or by an instance of
> Matrices.
I was just thinking while I was reading you : if he use generics, then
the generic package should define two matrix types, one for source
matrix and one for result matrix of some matrix functions (which
returns matrix with reversed dimensions). One Matrix(X,Y) and one
Matrix(Y,X)

On 2 f=E9v, 09:52, Jean-Pierre Rosen <ro...@adalog.fr> wrote:
> When designing a programming language, you have to stop at some point.
You have to stop to primitives (like rendezvous for tasking), things
which are a pain to simulate if not supported by the language (like
post/pre condition which will come with the next revision),
consistency (like my suggestion to make object declared with Constant
overloadable) and constructs with overall properties (like the loop
construct and its immutable variant instead of the goto-made-loop).

For every thing else (starting with personal-need-of-the-day), there
are generics and generic instantiations

(and a language is not a library, it's a paradigm)

On Wed, 3 Feb 2010 20:13:20 -0800 (PST), Hibou57 (Yannick Duch�ne) wrote:

> On 1 f�v, 03:11, "Peter C. Chapin" <pcc482...@gmail.com> wrote:
>> Apparently the component type of an array needs to be fully constrained (which
>> again makes sense) yet I don't know the size
> That's what generics are made for.
> But you will have a single matrix type of a fixed size.

Generics per design are incapable to express that the thing is variable.
You cannot make a generic variable-dimension matrix, a generic string, or
for that matter even a generic number:

generic
   Value : WHAT?
package Integer is
   ...
end Integer;

That does not work.

> I you don't want a fixed size type, the best I will think about, would
> be to create an abstract matrix type which will be implemented on a
> one dimensional array. You could use slice in the implementation.

This does not work either, I mean at the user interface level, because
slice does not have a type. In Ada types system you cannot express "the
subtype S is a vector of the type M". Therefore making it abstract types
you will loose most of the comfort built-in arrays offer.

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

On 4 f=E9v, 10:10, "Dmitry A. Kazakov" <mail...@dmitry-kazakov.de>
wrote:
> This does not work either, I mean at the user interface level, because
> slice does not have a type. In Ada types system you cannot express "the
> subtype S is a vector of the type M". Therefore making it abstract types
> you will loose most of the comfort built-in arrays offer.
I've was to post him some seed of a package (not tested, I sincerely
apologize for that) :

package Matrix is
   -- Provides a matrix type with an efficient implementation
   -- of a row swapping method. A similar method as the one
   -- provided here, may be used to get a matrix type with an
   -- efficient implementation of column swapping method.
   -- This package most likely need to be extended in order
   -- to be fully useful to something in real life.

   subtype Float_Type is Float;

   subtype One_Dimensional_Size_Type is Natural range 1 .. 2 ** 15;
   -- Either a number of rows or a number of columns.
   -- The maximum size is so that One_Dimensional_Size_Type ** 2 is
   -- most unlikely to commit an overflow and will mostly be a valid
   -- value on most of target machines.

   First_Index : constant :=3D 1;
   -- This may be changed to zero if needed, whenever the
   -- C indexing convention is preferred. This may not be
   -- changed to a value greater than one, otherwise
   -- Index_Type may potentially go out of Natural range.

   subtype Index_Type is Natural
      range
        (First_Index) ..
        (First_Index + (One_Dimensional_Size_Type'Last - 1));

   type Instance_Type (<>) is private;
   -- The type must be initialized at declaration :
   -- use New_Instance in that purpose (see below).

   -- The methods First_Row, Last_Row, First_Column,
   -- Last_Column and Item, are at least required to
   -- define pre- and post- condition of the New_Instance
   -- and Swap_Rows method.

   function First_Row
     (Instance : Instance_Type)
      return Index_Type;
   --|Ensures: Result =3D First_Index

   function Last_Row
     (Instance : Instance_Type)
      return Index_Type;

   function First_Column
     (Instance : Instance_Type)
      return Index_Type;
   --|Ensures: Result =3D First_Index

   function Last_Column
     (Instance : Instance_Type)
      return Index_Type;

   function Item
     (Instance : Instance_Type;
      Row      : Index_Type;
      Column   : Index_Type)
      return Float_Type;
   --|Requires: Row in First_Row (Instance) .. Last_Row (Instance);
   --|Requires:
   --|   Column in First_Column (Instance) ..
   --|   Last_Column (Instance);

   function New_Instance
      (Number_Of_Rows    : One_Dimensional_Size_Type;
       Number_Of_Columns : One_Dimensional_Size_Type)
       return Instance_Type;
   --|Ensures:
   --|   Last_Row (Instance) =3D
   --|   First_Index + Number_Of_Rows - 1;
   --|Ensures:
   --|   Last_Column (Instance) =3D
   --|   First_Index + Number_Of_Column - 1;
   --|Ensures:
   --|   for each Row in First_Row (Instance) .. Last_Row (Instance)
   --|   =3D> for each Column in First_Column (Instance) ..
   --|   Last_Column (Instance) =3D> Item (Row, Column) =3D 0.0;
   -- Note : all items are zero initialized.

   procedure Swap_Rows
      (Instance : in out Instance_Type;
       Row_1    : in Index_Type;
       Row_2    : in Index_Type);
   --|Requires: Row_1 in First_Row (Instance) .. Last_Row (Instance);
   --|Requires: Row_2 in First_Row (Instance) .. Last_Row (Instance);
   --|Ensures:
   --|   for each Column in First_Column (Instance) ..
   --|   Last_Column (Instance) =3D> (Item (Instance, Row_1, Column) =3D
   --|   old Item (Instance, Row_2, Column));
   --|Ensures:
   --|   for each Column in First_Column (Instance) ..
   --|   Last_Column (Instance) =3D> (Item (Instance, Row_2, Column) =3D
   --|   old Item (Instance, Row_1, Column));
   -- Note: Row_1 and Row_2 may be equal.

   --
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
   -- *** DO NOT CROSS *** DO NOT CROSS *** DO NOT CROSS *** DO NOT
CROSS
   --
-------------------------------------------------------------------

private

   pragma Assert
     ((First_Index + (One_Dimensional_Size_Type'Last ** 2)) <=3D
      (Natural'Last));
   -- Ensure computation on indexes will never go out of range.

   pragma Inline (First_Row);
   pragma Inline (Last_Row);
   pragma Inline (First_Column);
   pragma Inline (Last_Column);

   type Storage_Type is array (Natural range <>) of Float_Type;
   -- Flatenned array rows first columns next.
   -- Ex. (Row-1,Col-1),(Row-2,Col-1),(Row-1,Col-2),(Row-2,Col-2)
   --
   -- This way of flatening the matrix is the one which gives
   -- better performance to swap rows, as all items of a row
   -- appears to be consecutive, allowing slice access.
   --
   -- If good performance at swaping columns is targeting instead,
   -- then use the reverse representation (that is, like the
   -- Fortran convention), and update implementation of methods
   -- accordingly.
   --
   -- Note: the index type musn't be Index_Type,
   -- which is an index in one dimension, not in the overall storage.

   type Instance_Type
      (Number_Of_Rows    : One_Dimensional_Size_Type;
       Number_Of_Columns : One_Dimensional_Size_Type;
       Last_Data_Index   : Natural)
   is record
      Data : Storage_Type (First_Index .. Last_Data_Index);
   end record;

end Matrix;

package body Matrix is

   function First_Row
     (Instance : Instance_Type)
      return Index_Type
   is
   begin
      return First_Index;
   end First_Row;

   function Last_Row
     (Instance : Instance_Type)
      return Index_Type
   is
   begin
      return First_Index + (Instance.Number_Of_Rows - 1);
   end Last_Row;

   function First_Column
     (Instance : Instance_Type)
      return Index_Type
   is
   begin
      return First_Index;
   end First_Column;

   function Last_Column
     (Instance : Instance_Type)
      return Index_Type
   is
   begin
      return First_Index + (Instance.Number_Of_Columns - 1);
   end Last_Column;

   function Start_Of_Row
      (Instance : Instance_Type;
       Row     : Index_Type)
       return Natural
   -- Index in Data for the first item of the row
   -- whose index is Row. Used for accessing a matrix items
   -- and for row slice accesses.
   is
      Rows_Size : constant Natural :=3D Instance.Number_Of_Columns;
   begin
      return First_Index +  (Rows_Size * (Row - First_Index));
   end Start_Of_Row;

   pragma Inline (Start_Of_Row);

   function End_Of_Row
      (Instance : Instance_Type;
       Row     : Index_Type)
       return Natural
   -- Index in Data for the last item of the row
   -- whose index is Row. Used for row slice accesses.
   is
      Rows_Size : constant Natural :=3D Instance.Number_Of_Columns;
   begin
      return Start_Of_Row (Instance, Row) +  (Rows_Size - 1);
      -- This is an upper bound, not a limit, so
      -- Rows_Size - 1 is used. We add an offset starting
      -- the fist index of the row.
   end End_Of_Row;

   pragma Inline (End_Of_Row);

   function Item
     (Instance : Instance_Type;
      Row      : Index_Type;
      Column   : Index_Type)
      return Float_Type
   is
   begin
      Validity_Constraints :
      begin
         if Row not in First_Index .. Last_Row (Instance) then
            raise Program_Error;
         end if;
         if Column not in First_Index .. Last_Column (Instance) then
            raise Program_Error;
         end if;
      end Validity_Constraints;

      Method:
      declare
         Row_Start : constant Natural :=3D Start_Of_Row (Instance, Row);
         Column_Offset : constant Natural :=3D Column - First_Index;
         Data_Index : constant Natural :=3D Row_Start + Column_Offset;
      begin
         return Instance.Data (Data_Index);
      end Method;
   end Item; -- Procedure

   function New_Instance
      (Number_Of_Rows    : One_Dimensional_Size_Type;
       Number_Of_Columns : One_Dimensional_Size_Type)
       return Instance_Type
   is
      Data_Size : constant Natural :=3D
        (Number_Of_Rows * Number_Of_Columns);

      Last_Data_Index : constant Natural :=3D
        (First_Index + (Data_Size - 1));
   begin
      return
        (Number_Of_Rows    =3D> Number_Of_Rows,
         Number_Of_Columns =3D> Number_Of_Columns,
         Last_Data_Index   =3D> Last_Data_Index,
         Data              =3D> (others =3D> 0.0));
   end New_Instance; -- Function

   procedure Swap_Rows
      (Instance : in out Instance_Type;
       Row_1    : in Index_Type;
       Row_2    : in Index_Type)
   is
   begin
      Validity_Constraints :
      begin
         if Row_1 not in First_Index .. Last_Row (Instance) then
            raise Program_Error;
         end if;
         if Row_2 not in First_Index .. Last_Row (Instance) then
            raise Program_Error;
         end if;
      end Validity_Constraints;

      Special_Cases :
      begin
         if Row_1 =3D Row_2 then
            return;
         end if;
      end Special_Cases;

      Method :
      declare
         -- See comments in the package's private part.

         Rows_Size : constant Natural :=3D
           (Instance.Number_Of_Columns);

         Slice_Start_Of_Row_1 : constant Natural :=3D
            Start_Of_Row (Instance, Row_1);

         Slice_End_Of_Row_1 : constant Natural :=3D
            End_Of_Row (Instance, Row_1);

         Slice_Start_Of_Row_2 : constant Natural :=3D
            Start_Of_Row (Instance, Row_1);

         Slice_End_Of_Row_2 : constant Natural :=3D
            End_Of_Row (Instance, Row_2);

         Row_1_Slice : constant Storage_Type
         -- Backup of Row-1
           (Slice_Start_Of_Row_1 .. Slice_End_Of_Row_1) :=3D
               Instance.Data
                 (Slice_Start_Of_Row_1 ..
                  Slice_End_Of_Row_1);
      begin
         -- Copy Row-2 at the place of Row-1 (Row-1 was
         -- backed up).
         Instance.Data
           (Slice_Start_Of_Row_1 ..
            Slice_End_Of_Row_1)
            :=3D
               Instance.Data
                 (Slice_Start_Of_Row_2 ..
                  Slice_End_Of_Row_2);

         -- Copy backup of Row-1 at the place of
         -- Row-2.
         Instance.Data
           (Slice_Start_Of_Row_2 ..
            Slice_End_Of_Row_2)
            :=3D
               Row_1_Slice;
      end Method;

   end Swap_Rows; -- Procedure

end Matrix; -- Package

"Randy Brukardt" <randy@rrsoftware.com> writes:

> "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> wrote in message 
> news:178rg3rch8qdu$.13cgkywb09x3p.dlg@40tude.net...
>> On Sun, 31 Jan 2010 21:11:12 -0500, Peter C. Chapin wrote:
>>
>>> So I thought, "Perhaps A needs to be an array of arrays."
>>>
>>> type Matrix is array(Positive range <>) of WHAT_EXACTLY?
>>>
>>> Apparently the component type of an array needs to be fully constrained 
>>> (which
>>> again makes sense)
>>
>> Not really. Arrays should have been allowed to have discriminants. That 
>> was
>> missed in the language design.
>
> Early Ada 9x had discriminants for arrays. The capability got dropped when 
> lots of "nice-to-have" features got dropped from Ada 95. Such "array" types 
> have to be implemented as discriminant dependent records anyway; there is no 
> real hope of performance improvement from them, but a lot of implementation 
> complication. (I recall I was one of the stronger opponents of this feature, 
> mostly for implementation cost/benefit reasons.)
>
> So it is completely wrong to say that "this was missed in the language 
> design". The correct statement is that "this was explicitly rejected in the 
> language design".

It was missed in the design of Ada 83.  By the time of Ada 9X,
it was too late.  Sad.

If Ada 83 had had discriminated arrays, it would have been a
(slightly) simpler language.  And they're not hard to implement,
if you know about them ahead of time.  They may be hard to add
to an existing Ada 83 compiler that was designed based on
the rules of Ada 83.

That's true in general -- it's a lot easier to build an Ada 95 compiler,
than it is to build an Ada 83 compiler and then modify it to support
Ada 95.

Discriminated arrays would be more efficient than Ada's arrays.
We would have:

    type String (Length : Natural) is
        array (Positive range 1..Length) of Character;

This typically shrinks all unknown-length strings by 32 bits,
which is significant.  (Think about all the memory currently being
wasted to store zillions of copies of the number 1.)
And it makes bounds checking cheaper because the lower
bound is known at compile time.

And it would reduce the number of bugs -- all strings would
start at 1 (as opposed to the current sitation, where 99.99%
of all strings start at 1, and remaining 0.01% are bugs
waiting to happen).

And it would simplify the language (no need for the "range <>" syntax,
no need for the 'First and 'Last attributes of arrays, no need for
separate rules about array bounds and discriminants (which are
almost, but not quite, the same), etc).  Adding it to Ada 95 would NOT
simplify, because we'd have to keep all those now-useless features
for compatibility.

- Bob

On Feb 2, 8:52=A0am, Jean-Pierre Rosen <ro...@adalog.fr> wrote:
> Jerry a =E9crit :> I've never understood why Ada does not allow slicing i=
n
> > multidimensional arrays. What are the safety issues involved? And how
> > is it safe to force the programmer into ad hoc methods?
>
> One-dimensional slices are simple and efficient. Multidimensional slices
> are a can of worms.
>
> I guess you are thinking about rectangular slices. But why stop there? A
> slice comprising the main diagonal and the diagonal above and below can
> be very useful for some calculations. Or a slice which is the triangular
> part of the upper half...
>
> AFAICT, Fortran-99 does provide this - and the syntax is so complicated
> that nobody uses it. And implementation is also a nightmare.
>
> When designing a programming language, you have to stop at some point.
> The ratio (cost of implementation) / usefulness is a good measure for
> this. I think the ratio was simply to high for this feature.
>
> --
> ---------------------------------------------------------
> =A0 =A0 =A0 =A0 =A0 =A0J-P. Rosen (ro...@adalog.fr)
> Visit Adalog's web site athttp://www.adalog.fr


I've never felt the need for two dimensional (or higher
dimensional) slicing. It's partly a performance issue: if
you make the data storage matrix as small as possible (ie
make it exactly the same size as the matrix-transformation
you are doing) then you sometimes get a disastrous loss of
efficiency. If on the other hand you always make the data storage
matrix larger than necessary (so that your transformation will
always be performed on a sub-matrix of the data storage matrix),
then you always have the option of avoiding these efficiency
losses. Once you write the transformation routine so that it
operates on sub-matrices of the data storage matrix, then
you usually don't have to slice or copy a sub-matrix from
the data storage matrix in order to transform it.

Here are some Ada and Fortran examples of the
problem on various sized data storage arrays.
I used some bits and pieces of routines from
http://web.am.qub.ac.uk/users/j.parker/miscellany/

First example: we eigen-decompose an N x N =3D 2048 x 2048 matrix.
The data storage matrix is M x M =3D (1024+Padding) x (1024+Padding)
Here is the running time in seconds of an iterative jacobi
eigen-decomposition:

   2048x2048:  322 sec, gnat      (Padding=3D24)
   2048x2048: 1646 sec, gnat      (Padding=3D0)
   2048x2048: 1632 sec, gfortran  (Padding=3D0)
   2048x2048: 1492 sec, ifort     (Padding=3D0)

The observed 500% slowdown in the absence of padded arrays is
unacceptable, even if it is a rare event (occurring only on 2**p sized
data arrays). In fact it's not all that rare ... more comments on
that below.  (BTW, ifort is the INTEL fortran compiler, all
optimizations at Max. gfortran is the gcc variant.)

By the way, I never bothered to write a paddable Fortran routine,
but here are some more timings near a power-of-2 array bounds:

   1023x1023: 33.5 sec, gnat     (Padding=3D9)
   1023x1023: 42.4 sec, gfortran (Padding=3D0)
   1023x1023: 37.3 sec, ifort    (Padding=3D0)

   1022x1022: 33.5 sec, gnat     (Padding=3D10)
   1022x1022: 30.2 sec, gfortran (Padding=3D0)
   1022x1022: 28.3 sec, ifort    (Padding=3D0)

   1024x1024: 33.2 sec, gnat     (Padding=3D8)
   1024x1024: 96.0 sec, gnat     (Padding=3D0)
   1024x1024: 116  sec, gfortran (Padding=3D0)
   1024x1024: 43   sec, ifort    (Padding=3D0)

There is one puzzle here I don't have time solve.  Normally, a good
fortran will automatically pad the array for you ... I recall that
happening in the past. This time it seems to have slipped its fortran
mind. The ifort man pages:

   -O3 enables "Padding the size of certain power-of-two
   arrays to allow more  efficient cache use."

But I used -O3 and also -O3 -fast ... maybe did something wrong,
but important lesson: compiler optimization policies change
with time, and of course vary from compiler to compiler.
You can't rely on them or even, amazingly, man pages. It's much
better to write the program in a way that is insensitive to changing
optimization policies.

A web search of "cache thrashing" will reveal much depressing
detail on the subject. The efficiency problems discussed above
occur as our arrays become large and spill out of the L3 cache
(6 Meg in the present case).

Just to demonstrate that these problems show up on all sorts
of arrays, I did some runs in the 2000 to 3000 range, this time
using a householder decomposition scavenged from the Golub
singular value decomposition. Can still find plenty of 500%'s:

   2102x2102, 3.93 sec, gnat  (Padding =3D 0)
   2102x2102, 3.19 sec, gnat  (Padding =3D 8)

   2176x2176, 5.03 sec, gnat  (Padding =3D 0)
   2176x2176, 3.42 sec, gnat  (Padding =3D 8)

   2304x2304, 8.47 sec, gnat  (Padding =3D 0)
   2304x2304, 4.52 sec, gnat  (Padding =3D 8)

   2560x2560, 24.1 sec, gnat  (Padding =3D 0)
   2560x2560, 5.42 sec, gnat  (Padding =3D 8)

   3072x3072, 38.9 sec, gnat  (Padding =3D 0)
   3072x3072, 7.76 sec, gnat  (Padding =3D 8)

   3584x3584, 53.2 sec, gnat  (Padding =3D 0)
   3584x3584, 11.5 sec, gnat  (Padding =3D 8)

Finally, an example on 1-dim arrays, using fast fourier
transforms, FFT. The standard, and the most common FFT is
computed on a power-of-2 length data set: 0 .. 2**p-1. I
timed computation of these FFT's on arrays of length  2**p, and
I compared this with computation on arrays of
length  2**p + Padding, where Padding =3D 24.  The computation
on the padded arrays was faster. The ratio of running times is:

   p =3D 10  ratio =3D .93
   p =3D 11  ratio =3D .88
   p =3D 12  ratio =3D .76
   p =3D 13  ratio =3D .79
   p =3D 14  ratio =3D .76
   p =3D 15  ratio =3D .75
   p =3D 16  ratio =3D .77
   p =3D 17  ratio =3D .84
   p =3D 18  ratio =3D .75
   p =3D 19  ratio =3D .45
   p =3D 20  ratio =3D .63
   p =3D 21  ratio =3D .69
   p =3D 22  ratio =3D .69
   p =3D 22  ratio =3D .67
   p =3D 24  ratio =3D .62
   p =3D 25  ratio =3D .62

So the problem is more common here, and smaller. (These
efficiency losses I still consider unacceptable, especially
in a routine whose reason for existence is efficiency.)
The problem is still worse when you take FFTs of two
dimensional arrays.

There is another (and entirely independent) reason
I prefer routines that perform their transformations on
arbitrary sub-matrices (or on arbitrary diagonal blocks)
of the data storage matrix. After writing my 1st linear
algebra routine, I was very pleased with myself, but it
didn't really do what I wanted.  I realized I needed to
transform the diagonal sub-blocks of a large matrix,
and do it iteratively on arbitrarily sized diagonal
sub-blocks. It was a very simple matter to modify the code to
do this, and the result was so convenient that I've never
considered doing it otherwise.

Jonathan

Le Sun, 14 Feb 2010 01:42:13 +0100, jonathan <johnscpg@googlemail.com> a=
  =

=E9crit:
> First example: we eigen-decompose an N x N =3D 2048 x 2048 matrix.
> The data storage matrix is M x M =3D (1024+Padding) x (1024+Padding)
> Here is the running time in seconds of an iterative jacobi
> eigen-decomposition:
>
>    2048x2048:  322 sec, gnat      (Padding=3D24)
>    2048x2048: 1646 sec, gnat      (Padding=3D0)
>    2048x2048: 1632 sec, gfortran  (Padding=3D0)
>    2048x2048: 1492 sec, ifort     (Padding=3D0)
>
> The observed 500% slowdown in the absence of padded arrays is
> unacceptable, even if it is a rare event (occurring only on 2**p sized=

> data arrays). In fact it's not all that rare ... more comments on
> that below.  (BTW, ifort is the INTEL fortran compiler, all
> optimizations at Max. gfortran is the gcc variant.)
So this is mostly about representation clauses finally. Is that it ?

Do not know if you already know this document (as I remember I picked it=
  =

up from some one thread at comp.lang.ada), I've talked about on the othe=
r  =

fr.c.l.a :
http://research.scee.net/files/presentations/gcapaustralia09/Pitfalls_of=
_Object_Oriented_Programming_GCAP_09.pdf
I had pointed about frames  #17, #18, #19 et #20, which contains good  =

source of inspiration. Hope this could help you to figure a path.

You've posted a long list of tests-bench and observations. I did not  =

looked at every thing, but hope I will have a more closer look at it lat=
er.

-- =

No-no, this isn't an oops ...or I hope (TM) - Don't blame me... I'm just=
  =

not lucky

On Feb 14, 1:54=A0am, Hibou57 (Yannick Duch=EAne)
<yannick_duch...@yahoo.fr> wrote:
> Le Sun, 14 Feb 2010 01:42:13 +0100, jonathan <johns...@googlemail.com> a =
=A0
> =E9crit:

  First example: we eigen-decompose an N x N =3D 2048 x 2048 matrix.
  The data storage matrix is M x M =3D (1024+Padding) x (1024+Padding)

should be of course:

  First example: we eigen-decompose an N x N =3D 2048 x 2048 matrix.
  The data storage matrix is M x M =3D (2048+Padding) x (2048+Padding)

> Do not know if you already know this document (as I remember I picked it =
=A0
> up from some one thread at comp.lang.ada), I've talked about on the other=
 =A0
> fr.c.l.a :http://research.scee.net/files/presentations/gcapaustralia09/Pi=
tfalls...
> I had pointed about frames =A0#17, #18, #19 et #20, which contains good =
=A0
> source of inspiration. Hope this could help you to figure a path.

Yes, I remembered this, probably from an old post of yours.  I wanted
to cite it when
when I posted earlier, but I could not find the site.  This is not
something you
forget quickly (frames 17 and 18):

   1980: RAM latency ~ 1 cycle
   2009: RAM latency ~ 400+ cycles

It's the heart of the matter, and it is just getting worse.  Helps
convince me
anyway that I did not waste time on an unimportant matter! In
numerical linear
algebra the usual solution is to restructure matrices as a collection
of
blocks.  That has a few of its own problems though. Minor footnote: I
did
some tests on Intel's new nehalem CPU's.  Vastly improved performance
on these
multi-megabyte arrays.  Problem not cured though. Don't know enough to
say more about it.
Thanks for the reminder.

Jonathan

On Tue, 02 Feb 2010 09:52:31 +0100, Jean-Pierre Rosen
<rosen@adalog.fr> wrote:

> Jerry a �crit :
> > I've never understood why Ada does not allow slicing in
> > multidimensional arrays. What are the safety issues involved? And how
> > is it safe to force the programmer into ad hoc methods?
> > 
> One-dimensional slices are simple and efficient. Multidimensional slices
> are a can of worms.
> 
> I guess you are thinking about rectangular slices. But why stop there? A
> slice comprising the main diagonal and the diagonal above and below can
> be very useful for some calculations. Or a slice which is the triangular
> part of the upper half...
> 
> AFAICT, Fortran-99 does provide this - and the syntax is so complicated
> that nobody uses it. And implementation is also a nightmare.
> 
Fortran 90 (and later) has rectangular but optionally non-unit-stride
slices; X(1:5:4,2:6:2) is X(1,2) X(1,4) X(1,6) X(5,2) X(5,4) X(5,6).
(Fortran arrays are column-major). (And although it treats string --
fixed length only -- as a different type than array of character, you
can use corresponding substrings of elements of an array of strings in
a way that is effectively also rectangular.)

It also has 'vector subscripting;: X(Y) accesses X(Y(1)) X(Y(2)) ...
X(Y(N)) -- but this cannot be used as modifiable actual argument or
POINTER target, making it not really a firstclass object/variable.

A major new 'coarray' feature, which essentially allows distribution
across parallel processing systems, vaguely like an in-language OpenMP
(although I hear the standardese must be rather contorted to avoid
impermissibly specifying implementation) was proposed for what was
scheduled to be F 08, but proved contentious enough that AFAIK it
still hasn't been finalized.

*PL/I* has the X DEFINED Y (iSUB) syntax which allows things like a
diagonal (but AFAICS not more than one at a time).

> When designing a programming language, you have to stop at some point.
> The ratio (cost of implementation) / usefulness is a good measure for
> this. I think the ratio was simply to high for this feature.

The features in F90 at least in this area weren't too much of a
problem, at least judging from the reports of visibly intelligent and
apparently informed people in c.l.f. Although those implementors had
the advantage that F90 was originally scheduled for I believe 87, and
even before that there had been experience with nonstandard but fairly
widespread HPF "High Performance Fortran" extensions.

In contrast F03, with adds features mostly in the areas of OO and
'parameterized' types somewhat like Ada discriminated ones, has taken
longer although most vendors are now reportedly getting close.

jonathan <johnscpg@googlemail.com> writes:

<snip>
> It's the heart of the matter, and it is just getting worse.  Helps
> convince me
> anyway that I did not waste time on an unimportant matter! In
> numerical linear
> algebra the usual solution is to restructure matrices as a collection
> of
> blocks.  That has a few of its own problems though. Minor footnote: I
> did

Another approach is to keep the matrix as a single big matrix, but
reformulate your loops as recursive procedure which divides the matrix
a number of times and then loops at the lowest level. The increased
efficiency from tha cache will typically more than outweigh the
overhead from the recursion.

The following is a java version of the idea applied to transposing a
lagre matrix. (I have Ada and C versions as well) The recursive
version is typically about ten times faster than the pure nested loop
version. It is also interesting to see the speed of tha java version
relative to a C or Ada version. Try it and see...


import java.util.*;

class transpose
{
    private static final int u1 = 6000;
    private static final int u2 = 6000;
    private static int[][] a = new int[u1][u2];
    private static int[][] b = new int[u1][u2];

    private static void init(int[][]x)
    {
        int i = 0;
        int j = 0;
        for (i = 0; i < x.length; i++) {
            for (j = 0; j < x[i].length; j++) {
                x[i][j] = i; 
            }
        }
    }

    private static void loop_transpose(int[][]a, int[][]b, int ll1, int ul1, int ll2, int ul2)
    {
        int i = 0;
        int j = 0;
        for (i = ll1; i <= ul1; i++) {
            for (j = ll2; j <= ul2; j++) {
                b[j][i] = a[i][j]; 
            }
        }
    }

    private static void recursive_transpose(int[][]a, int[][]b, int ll1, int ul1, int ll2, int ul2, int cutoff) 
    {
        int split = 0;
        if (ul1 - ll1 > cutoff || ul2 - ll2 > cutoff) {
            if (ul1 - ll1 > ul2 - ll2 ) {
                split = (ul1 - ll1)/2;
                recursive_transpose(a,b,ll1,ll1+split,ll2,ul2,cutoff);
                recursive_transpose(a,b,ll1+1+split,ul1,ll2,ul2,cutoff);
            }else{
                split = (ul2 - ll2)/2;
                recursive_transpose(a,b,ll1,ul1,ll2,ll2+split,cutoff);
                recursive_transpose(a,b,ll1,ul1,ll2+1+split,ul2,cutoff);
            }
        }else {
            loop_transpose(a,b,ll1,ul1,ll2,ul2);
        }
    }

    public static void main(String[] args)
    {
        init(a); init(b);

        long start = (new Date()).getTime();
        loop_transpose(a,b,0,u1-1,0,u2-1);

        long stop = (new Date()).getTime();
        System.out.println( "loop " + (stop-start)); 

        start = (new Date()).getTime();
        recursive_transpose(a,b,0,u1-1,0,u2-1,16);

        stop = (new Date()).getTime();
        System.out.println( "recursive " + (stop-start)); 
    }
}
 
-- 
   C++: The power, elegance and simplicity of a hand grenade.