Re: on using array index, vectorized operation

[Phil Clayton <phil.clayton@lineone.net>]

On Mar 29, 5:55 am, "Nasser M. Abbasi" <n...@12000.org> wrote:
> I rewrote the procedure using this package (I am using float, so
> it worked for me), and was really surprised that now I able
> to remove the loop and use vectorized statement:
>
> -------------
>        subtype r is Natural range o.u'First + 1 .. o.u'Last - 1;
>        subtype r_plus_1 is Natural range r'First + 1 .. r'Last + 1;
>        subtype r_minus_1 is Integer range r'First - 1 .. r'Last - 1;
>
>     begin
>
>        u (r) := u (r) -
>                 (a / 2.0) * (u (r_plus_1) - u (r_minus_1)) +
>                 (a ** 2) / 2.0 *
>                 (u (r_minus_1) - 2.0 * u (r) + u(r_plus_1));
>
> -------------
>
> It is now almost the same form as in Fortran and Matlab !
>
> Only extra thing is the need to define those subtype's for
> the ranges. (r, r_plus_1, r_minus_1) .
>
> In Fortran, I would have written the above as
>
> -----------------------------------------
>        u (i) := u (i) -
>                 (a / 2.0) * (u (i+1) - u (i-1)) +
>                 (a ** 2) / 2.0 *
>                 (u (i-1) - 2.0 * u (i) + u(i+1) )
> ----------------------------------------------
>
> Where 'i' is an array of the index range I am interested in.

It's worth observing that you're using vectorized subscripting, not
just vectorized arithmetic.  It should be possible to write a function
that does the vectorized subscript, e.g. Sub(A, X) = B such that B(i)
= A(X(i)).  Unfortunately, if you use real arrays for the i+1 in your
example, you'll have real subscripts.  That may be acceptable when
modelling e.g. in Matlab but not in an implementation.  So you may
need + and - on arrays of your index type.  It all depends what you
intend to use your code for.

Phil

P.S. I still think Ada should allow us to write something like

  u (2 .. 5) := (for i in 2 .. 5 =>
                  u (i) -
                  (a / 2.0) * (u (i+1) - u (i-1)) +
                  (a ** 2) / 2.0 *
                  (u (i-1) - 2.0 * u (i) + u(i+1)));

(which avoids the temporary variable that the equivalent loop
requires).  This enables on-the-fly vectorization by promoting scalar
operations inline in expressions.  For example, if you didn't have a
version of some function F vectorized in its second argument, you
could write

  (for I in B'Range => F(A, B(I), C))