questions on using the array component iterator for 2012

questions on using the array component iterator for 2012

[Nasser M. Abbasi]

reference

http://www.ada-auth.org/standards/12rm/html/RM-5-5-2.html

The  array component iterator seems nice, but I was not sure
if it can be used access elements in the matrix that are
indexed off the current iterator.

For example, assume I write

---------------------------------------------
with Ada.Numerics.Real_Arrays;  use Ada.Numerics.Real_Arrays;
procedure foo is
    A : constant Real_Matrix :=
          ( ( 1.0,  2.0, 3.0 ) ,
            ( 3.0,  5.0, 6.0 ) ,
            ( 7.0,  8.0, 9.0 ) );
begin
    FOR e of A  LOOP
        e := e + 1 ;
    END
end Matrix_Product;
---------------------------------------------

So the above adds '1' to each element in the matrix.
(order of iterations is defined in the RM, say I used
FORTRAN order).

But suppose I want to implement say a Jacobi relaxation
iteration to solve a 2D Laplace pde on that grid, which
is defined to update each element in the grid using

A(i) := (1/4) (  A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1)

where here 'i' is row index, and 'j' is column index.

http://en.wikipedia.org/wiki/Relaxation_%28iterative_method%29

(iteration has to start from the second row and second column,
to avoid falling off the edge of the grid ofcourse, but this
is not important now)

i.e the above just replaces the current element in the matrix by
the average of the 4 adjacent elements in the little grid around
the current element.

My question is:

using 'e', this new iterator, can one somehow implement the above,
by somehow referencing these 4 adjacent elements in the grid?
(above, left, right, below)?

My point is that, if the iterator can be used only to access
the current element, then this might not be too useful.

For example, in the above example toy example, I do not
even need an iterator, I can use direct operation in Fortran or
Matlab say and just type

           A = 1.0 + A

(I do not know if one can do this in Ada btw?) May be
one needs to just define a "+" function for this. (I am not
too experienced in Ada, will try it).

Thanks
--Nasser

Re: questions on using the array component iterator for 2012

[Nasser M. Abbasi]

On 6/23/2012 7:02 AM, Nasser M. Abbasi wrote:

> ---------------------------------------------
> with Ada.Numerics.Real_Arrays;  use Ada.Numerics.Real_Arrays;
> procedure foo is
>      A : constant Real_Matrix :=
>            ( ( 1.0,  2.0, 3.0 ) ,
>              ( 3.0,  5.0, 6.0 ) ,
>              ( 7.0,  8.0, 9.0 ) );
> begin
>      FOR e of A  LOOP
>          e := e + 1 ;
>      END
> end Matrix_Product;
> ---------------------------------------------
>

opps, need to remove the 'constant' word above, so that A can be updated
ofocurse.

--Nasser

Re: questions on using the array component iterator for 2012

[Adam Beneschan]

On Saturday, June 23, 2012 5:02:18 AM UTC-7, Nasser M. Abbasi wrote:
> reference

> 
> But suppose I want to implement say a Jacobi relaxation
> iteration to solve a 2D Laplace pde on that grid, which
> is defined to update each element in the grid using
> 
> A(i) := (1/4) (  A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1)
> 
> where here 'i' is row index, and 'j' is column index.
> 
> http://en.wikipedia.org/wiki/Relaxation_%28iterative_method%29

That wouldn't even work, I think.  My understanding of matrix math is limited, but the Wikipedia article talks about repeatedly performing a *matrix* assignment "phi := phi*" where the elements of phi* are defined in terms of the elements of phi.  That is, you have to compute the entire matrix phi* first, before reassigning back to the original matrix.  The code you've written above will modify elements of A before they need to be used.  You'd need to do something like

New_A(i,j) := (1/4) (  A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1)

(you left the ",j" off the index of the left-hand side)

and then A := New_A after all elements are computed.

At least that's how it looks to me.  I don't really understand the math here, and maybe it's possible that the end results, if the loop is repeated enough times, would be the same even if you do modify elements of A before they're used.  But in any event, (1) an Ada array iterator isn't going to work if more than one array is involved, and (2) I think array iterators really are intended only for when you're focusing on each element one at a time.  If an operation on an array element requires that you look at other elements of the same array, then an array iterator isn't an appropriate construct to use.  You'd want to make the indexes (in your case, the row and column numbers) explicit so that readers will understand what you're doing.  

                              -- Adam

Re: questions on using the array component iterator for 2012

[Nasser M. Abbasi]

On 6/25/2012 1:19 PM, Adam Beneschan wrote:
> On Saturday, June 23, 2012 5:02:18 AM UTC-7, Nasser M. Abbasi wrote:

> That wouldn't even work, I think.
...
>
> New_A(i,j) := (1/4) (  A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1)
>

For the LHS, there is  a jacobi iteration, which needs a
temporary matrix, like you have above, and there is gauss-siedel
iteration, which uses the same (in place) matrix. There are
many other types of iterations (called relaxations methods such
as sor, gauss-siedel with red/black variation, sor with chebyshev,
and may be more).

But that is all not really important.

What I was asking for, is using the new iterator, if one can
use it to access elements outside its "sphere of influence",
ie. to access 1-off it, or 2-off of it, up and down, etc. i.e.
convert it to an index locally so as to visits elements other
than it somehow.

But after thinking about it, I do not think this is
practical to do. If I wanted to do the above type iteration,
I'll just normal loops and use normal indices, not the
new iterator.

> (you left the ",j" off the index of the left-hand side)
>

yes, sure. I was writing things on the fly.

> and then A := New_A after all elements are computed.
>

Yes. For Jacobi. I've written all these algorithms many times
before in Matlab and in Mathematica. If you interested in the code,
Check my Mathematica applets here

http://12000.org/my_notes/mma_demos/index.htm

#25 in the table implements many solvers including Jacboi.
You can run in it the browser. (need plugin, free of charge,
see link there if interested). Firefox, etc.. but no plugin for
Linux yet, but free player yes, can be downloaded from Wolfram
research.  

regards,
--Nasser

Re: questions on using the array component iterator for 2012

[Shark8]

On Saturday, June 23, 2012 7:02:18 AM UTC-5, Nasser M. Abbasi wrote:
> For example, assume I write
>     FOR e of A  LOOP
>         e := e + 1 ;
>     END LOOP
> [...]
> 
> But suppose I want to implement say a Jacobi relaxation
> iteration to solve a 2D Laplace pde on that grid, which
> is defined to update each element in the grid using
> 
> A(i,j) := (1/4) (  A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1)
> 
> where here 'i' is row index, and 'j' is column index.

Well the problem is that the "FOR E IN ..." loop is conceptually different from a "FOR Index IN Array_Variable'RANGE" loop. The difference is that in the former you don't care what the index/indices are, you just want to apply the same procedure to each element.

Given that needing the index is sometimes wanted, the designers may have added a 'Index(X) attribute, I don't really know because I've not had the inclination to explore the new loop in such depth.