We all know what

a(:,:)=b(:,:)+c(:,:)

will do

but what about

real ul(8,3), uu(10,20,3)

integer itabsrt(8,3)

ul(:,:)=uu(itabsrt(:,1),itabsrt(:,2),itabsrt(:,3),:)       ! legal?

The intention is

do i=1,8
do j=1,3
    ul(i,j)=uu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j)
enddo
enddo

It's probably not allowed but I can't think of a specific reason.

OR how about

ul(:,:)=uu(itabsrt(:,:),itabsrt(:,:),itabsrt(:,:),2)          ! legal?

It seems to me this is just lke the a=b+c example,

OR, with

real ur(10,10)

ur(:,:)=uu(itabsrt(:,1),:)                                       !
legal?

Basically array operations can get very complicated especially when
the array indices are themselves found from array references. None of
the books I've read really spell out formal rules for arbitrarily
complex array operations.

On Oct 5, 12:44=A0pm, The Star King <j...@npl.co.uk> wrote:
> We all know what
>
> a(:,:)=3Db(:,:)+c(:,:)
>
> will do
>
> but what about
>
> real ul(8,3), uu(10,20,3)
>
> integer itabsrt(8,3)
>
> ul(:,:)=3Duu(itabsrt(:,1),itabsrt(:,2),itabsrt(:,3),:) =A0 =A0 =A0 ! lega=
l?

No, this is not legal, because you're associating a rank-4 array on
the r.h.s. to a rank-2 array in the l.h.s.! Also, you defined your
array 'uu' as rank-3, but it's being called as rank-4! There must be a
typo there somewhere, so it's important that you present your problem
more clearly. I'll assume the last colon shouldn't be there:

ul(:,:)=3Duu(itabsrt(:,1),itabsrt(:,2),itabsrt(:,3))

It's still not legal because this associates a rank-3 on the r.h.s. to
a rank-2 in the l.h.s.!

> The intention is
>
> do i=3D1,8
> do j=3D1,3
> =A0 =A0 ul(i,j)=3Duu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j)
> enddo
> enddo
>
> It's probably not allowed but I can't think of a specific reason.

Now I'm confused with the rank of 'uu'... still, assuming 'uu' is a
rank-4 and that you forgot to declare one of the dimensions, all of
what I said above is still valid. You cannot associate a rank-2 array
on the left with a rank-4 array on the right! So your nested DO loop
is also not legal!

> OR how about
>
> ul(:,:)=3Duu(itabsrt(:,:),itabsrt(:,:),itabsrt(:,:),2) =A0 =A0 =A0 =A0 =
=A0! legal?

Even worse! You're associating a rank-2 with a rank-6!

> It seems to me this is just lke the a=3Db+c example,

It is not! I'll explain bellow.

> OR, with
>
> real ur(10,10)
>
> ur(:,:)=3Duu(itabsrt(:,1),:) =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 !
> legal?

Finally! Rank-2 on the left and rank-2 on the right. This is a good
start, but not enough. It is only a legal association if the
components of ITABSRT are each is an integer comprised between 1 and
10 (the first dimension of 'uu')!

> Basically array operations can get very complicated especially when
> the array indices are themselves found from array references. None of
> the books I've read really spell out formal rules for arbitrarily
> complex array operations.

No, these are not complex. I just think you didn't fully understand or
ever came across the concept of vector subscripts. ITABSRT is just an
integer array that is supposed to contain valid integer indices of
'ur' or 'ul'.

Consider this simple example:

integer :: v(4)
real    :: a(9)
real    :: b(4)

v =3D (/ 4, 3, 1, 7 /)
b =3D a(v)

This is the same as doing:

b(:) =3D a(v(:))

or:

b =3D (/ a(4), a(3), a(1), a(7) /)

See what's going on? The colon inside v(:) simply says that all
components of 'v' will be used as indices of 'a'; these indices
(integer components) ought to be less or equal than the dimension of
'a', obviously! For example:

v =3D (/ 4, 3, 1, 10 /)

would give a segmentation fault, because 10 > 9 (and a(10) is out of
array bounds). On the other hand, the colon inside b(:) says that 'b'
will collect all elements of 'a' whose indices are given by 'v', and
in that order. So, another requirement for a valid association is that
'v' and 'b' have the same dimension.

AFAIK, only rank-1 array (vector) subscripts are possible in Fortran.
You just need to be very careful how you manipulate them. If you
follow my simple example above, you'll see that everything else is a
very simple generalization.

And try to present your problem more clearly and carefully, because
there is clearly a lot of confusion in your examples, regarding array
dimensions, and in the end I didn't understand what is your objective.

Best,

-- helvio

On Oct 5, 1:45=A0pm, helvio <helvio.vairin...@googlemail.com> wrote:
> On Oct 5, 12:44=A0pm, The Star King <j...@npl.co.uk> wrote:
>
> > We all know what
>
> > a(:,:)=3Db(:,:)+c(:,:)
>
> > will do
>
> > but what about
>
> > real ul(8,3), uu(10,20,3)
>
> > integer itabsrt(8,3)
>
> > ul(:,:)=3Duu(itabsrt(:,1),itabsrt(:,2),itabsrt(:,3),:) =A0 =A0 =A0 ! le=
gal?
>
> No, this is not legal, because you're associating a rank-4 array on
> the r.h.s. to a rank-2 array in the l.h.s.! Also, you defined your
> array 'uu' as rank-3, but it's being called as rank-4! There must be a
> typo there somewhere, so it's important that you present your problem
> more clearly. I'll assume the last colon shouldn't be there:
>
> ul(:,:)=3Duu(itabsrt(:,1),itabsrt(:,2),itabsrt(:,3))
>
> It's still not legal because this associates a rank-3 on the r.h.s. to
> a rank-2 in the l.h.s.!
>
> > The intention is
>
> > do i=3D1,8
> > do j=3D1,3
> > =A0 =A0 ul(i,j)=3Duu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j)
> > enddo
> > enddo
>
> > It's probably not allowed but I can't think of a specific reason.
>
> Now I'm confused with the rank of 'uu'... still, assuming 'uu' is a
> rank-4 and that you forgot to declare one of the dimensions, all of
> what I said above is still valid. You cannot associate a rank-2 array
> on the left with a rank-4 array on the right! So your nested DO loop
> is also not legal!
>
> > OR how about
>
> > ul(:,:)=3Duu(itabsrt(:,:),itabsrt(:,:),itabsrt(:,:),2) =A0 =A0 =A0 =A0 =
=A0! legal?
>
> Even worse! You're associating a rank-2 with a rank-6!
>
> > It seems to me this is just lke the a=3Db+c example,
>
> It is not! I'll explain bellow.
>
> > OR, with
>
> > real ur(10,10)
>
> > ur(:,:)=3Duu(itabsrt(:,1),:) =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 !
> > legal?
>
> Finally! Rank-2 on the left and rank-2 on the right. This is a good
> start, but not enough. It is only a legal association if the
> components of ITABSRT are each is an integer comprised between 1 and
> 10 (the first dimension of 'uu')!
>
> > Basically array operations can get very complicated especially when
> > the array indices are themselves found from array references. None of
> > the books I've read really spell out formal rules for arbitrarily
> > complex array operations.
>
> No, these are not complex. I just think you didn't fully understand or
> ever came across the concept of vector subscripts. ITABSRT is just an
> integer array that is supposed to contain valid integer indices of
> 'ur' or 'ul'.
>
> Consider this simple example:
>
> integer :: v(4)
> real =A0 =A0:: a(9)
> real =A0 =A0:: b(4)
>
> v =3D (/ 4, 3, 1, 7 /)
> b =3D a(v)
>
> This is the same as doing:
>
> b(:) =3D a(v(:))
>
> or:
>
> b =3D (/ a(4), a(3), a(1), a(7) /)
>
> See what's going on? The colon inside v(:) simply says that all
> components of 'v' will be used as indices of 'a'; these indices
> (integer components) ought to be less or equal than the dimension of
> 'a', obviously! For example:
>
> v =3D (/ 4, 3, 1, 10 /)
>
> would give a segmentation fault, because 10 > 9 (and a(10) is out of
> array bounds). On the other hand, the colon inside b(:) says that 'b'
> will collect all elements of 'a' whose indices are given by 'v', and
> in that order. So, another requirement for a valid association is that
> 'v' and 'b' have the same dimension.
>
> AFAIK, only rank-1 array (vector) subscripts are possible in Fortran.
> You just need to be very careful how you manipulate them. If you
> follow my simple example above, you'll see that everything else is a
> very simple generalization.
>
> And try to present your problem more clearly and carefully, because
> there is clearly a lot of confusion in your examples, regarding array
> dimensions, and in the end I didn't understand what is your objective.
>
> Best,
>
> -- helvio

Sorry! I did make a mistake. The declarations should be

real ul(8,3), uu(10,20,3,4)
integer itabsrt(8,3)

that is, uu should be rank 4.

I'm sure my do loops, at least, are legal!

do i=3D1,8
do j=3D1,3
    ul(i,j)=3Duu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j)
enddo
enddo

itabsrt(i,1) is just a scalar integer. I can use it to index uu.

uu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j) is a scalar real. I am
setting ul(i,j) (a scalar) to this number

Your remarks seem to indicates the ranks of array on left and right
must be the same. I don't think this is necessaryilu the case, eg

a(:)=3Db(:,3)

the array on the left is rank 1 and the right is rank 2, but because I
am taking an appropriate slice I can do it. (if this is illegal, I
REALLY am confused!). That's what I was attempting to do in my
examples.

The Star King wrote:
....

> Your remarks seem to indicates the ranks of array on left and right
> must be the same. I don't think this is necessarily the case, eg
> 
> a(:)=b(:,3)
> 
> the array on the left is rank 1 and the right is rank 2, but because I
> am taking an appropriate slice I can do it. (if this is illegal, I
> REALLY am confused!). That's what I was attempting to do in my
> examples.

The problem in your thinking is that the array on the RHS _IS_ rank one 
because the rhs is evaluated and then assigned to the lhs.  The array b 
doesn't enter into it in the sense other than you must write a proper 
indexing expression for it which is what helvio was saying that in your 
original posting there was a mismatch there.

--

On Oct 5, 3:32=A0pm, dpb <n...@non.net> wrote:
> The Star King wrote:
>
> ...
>
> > Your remarks seem to indicates the ranks of array on left and right
> > must be the same. I don't think this is necessarily the case, eg
>
> > a(:)=3Db(:,3)
>
> > the array on the left is rank 1 and the right is rank 2, but because I
> > am taking an appropriate slice I can do it. (if this is illegal, I
> > REALLY am confused!). That's what I was attempting to do in my
> > examples.
>
> The problem in your thinking is that the array on the RHS _IS_ rank one
> because the rhs is evaluated and then assigned to the lhs. =A0The array b
> doesn't enter into it in the sense other than you must write a proper
> indexing expression for it which is what helvio was saying that in your
> original posting there was a mismatch there.
>
> --

Thank you both, I am starting to understand now. So this

ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)

would be legal? I assume the mapping in colons is left to right. So
the left hand colon in ul(:,:) maps to the itabsrt(:,1) colon? So the
equivalent would be

do i=3D1,8
  do j=3D1,3
     ul(i,j)=3Duu(3,itabsrt(i,2),4,j)
  enddo
enddo

And I can also do this?

ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)+itabsrt(:,:)

giving

do i=3D1,8
  do j=3D1,3
     ul(i,j)=3Duu(3,itabsrt(i,2),4,j)+itabsrt(i,j)
  enddo
enddo

The Star King wrote:
....

> ul(1:8,:)=uu(3,itabsrt(1:8,2),4,:)
 > would be legal? ...

I think so; at least I didn't catch a problem...

> ... I assume the mapping in colons is left to right. So
> the left hand colon in ul(:,:) maps to the itabsrt(:,1) colon? ...

Again, _NO_!!!

The rhs is evaluated _COMPLETELY_in_ignorance_ of the LHS and that 
result is then _ASSIGNED_ to the LHS.

That those two entities must be conformant should be obvious, but the 
one "knows nuthink'" (to para-quote Sgt Schultz) of the other.  This 
isn't anything specific about or to array slices; it's the fundamental 
way Fortran works in any assignment.

--

On Oct 5, 4:55=A0pm, dpb <n...@non.net> wrote:
> The Star King wrote:
>
I think you need to refer to a Fortran text, as a minumimum
http://en.wikipedia.org/wiki/Fortran_language_features#Array_subobjects_.28=
sections.29
to find out about rank and conformance. Note, in particular, that
given an array a(10, 10) that

a(1, 1)     is a scalar (rank zero)
a(1:1, 1)   is an array section (rank one)
a(1:1, 1:1) is an array section (rank two)

Regards,

Mike Metcalf

On Oct 5, 4:24=A0pm, m_b_metcalf <michaelmetc...@compuserve.com> wrote:
> On Oct 5, 4:55=A0pm, dpb <n...@non.net> wrote:> The Star King wrote:
>
> I think you need to refer to a Fortran text, as a minumimumhttp://en.wiki=
pedia.org/wiki/Fortran_language_features#Array_subobjec...
> to find out about rank and conformance. Note, in particular, that
> given an array a(10, 10) that
>
> a(1, 1) =A0 =A0 is a scalar (rank zero)
> a(1:1, 1) =A0 is an array section (rank one)
> a(1:1, 1:1) is an array section (rank two)
>
> Regards,
>
> Mike Metcalf

Thanks, Mike, I've actually got (and read!) your book Fortran 95/2003
explained. I also attended a f90 lecture series you gave in CERN many
years ago...!

(so I'm probably pretty dumb)

On the specific point, may I ask, would the statement,

ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)

have the same effect as,

do i=3D1,8
  do j=3D1,3
     ul(i,j)=3Duu(3,itabsrt(i,2),4,j)
  enddo
enddo

?

In response to dpb:

> Again, _NO_!!!

> The rhs is evaluated _COMPLETELY_in_ignorance_ of the LHS and that
> result is then _ASSIGNED_ to the LHS.

This may be conceptually how one should think of it. But what the
compiler may actually do is convert matrix operations into do loops
and then compile. This is just how g95 works, I asked Andy Vaught
about it. That is, when confronted with

ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)

g95 will not create temporary arrays. (it will create temporaries if
the left side variable appears also on the right hand side)

On Oct 5, 5:37=A0pm, The Star King <j...@npl.co.uk> wrote:
> On Oct 5, 4:24=A0pm, m_b_metcalf <michaelmetc...@compuserve.com> wrote:
>
> > On Oct 5, 4:55=A0pm, dpb <n...@non.net> wrote:> The Star King wrote:
>
> > I think you need to refer to a Fortran text, as a minumimumhttp://en.wi=
kipedia.org/wiki/Fortran_language_features#Array_subobjec...
> > to find out about rank and conformance. Note, in particular, that
> > given an array a(10, 10) that
>
> > a(1, 1) =A0 =A0 is a scalar (rank zero)
> > a(1:1, 1) =A0 is an array section (rank one)
> > a(1:1, 1:1) is an array section (rank two)
>
> > Regards,
>
> > Mike Metcalf
>
> Thanks, Mike, I've actually got (and read!) your book Fortran 95/2003
> explained. I also attended a f90 lecture series you gave in CERN many
> years ago...!
>
> (so I'm probably pretty dumb)
>
> On the specific point, may I ask, would the statement,
>
> ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)
>
> have the same effect as,
>
> do i=3D1,8
> =A0 do j=3D1,3
> =A0 =A0 =A0ul(i,j)=3Duu(3,itabsrt(i,2),4,j)
> =A0 enddo
> enddo
>
> ?
>
> In response to dpb:
>
> > Again, _NO_!!!
> > The rhs is evaluated _COMPLETELY_in_ignorance_ of the LHS and that
> > result is then _ASSIGNED_ to the LHS.
>
> This may be conceptually how one should think of it. But what the
> compiler may actually do is convert matrix operations into do loops
> and then compile. This is just how g95 works, I asked Andy Vaught
> about it. That is, when confronted with
>
> ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)
>
> g95 will not create temporary arrays. (it will create temporaries if
> the left side variable appears also on the right hand side)

No, because the expressions don't match the array declaration. This
works:

real :: ul(8,3), uu(10,20,3, 4) =3D 1
integer ::itabsrt(8,3)
itabsrt(:8, 2) =3D (/ (i, i =3D 1, 8) /)
ul(1:8,:)=3Duu(3,itabsrt(1:8,2),:,4)
print *, ul
do i=3D1,8
  do j=3D1,3
     ul(i,j)=3Duu(3,itabsrt(i,2),j,4)
  enddo
enddo
print*, ul
end

On your point about how a compiler handles it in practice, that is of
no relevance since the standard makes no statement about how any
program should be compiled and executed, other than it be done
correctly.

Nice to hear from you about distant times. Since you have MR&C, please
(re)read Section 6.13.

Hope that helps,

Mike Metcalf

The Star King <jfb@npl.co.uk> wrote:

> > The rhs is evaluated _COMPLETELY_in_ignorance_ of the LHS and that
> > result is then _ASSIGNED_ to the LHS.
> 
> This may be conceptually how one should think of it. But what the
> compiler may actually do is convert matrix operations into do loops
> and then compile. This is just how g95 works, I asked Andy Vaught
> about it. That is, when confronted with
> 
> ul(1:8,:)=uu(3,itabsrt(1:8,2),4,:)
> 
> g95 will not create temporary arrays. (it will create temporaries if
> the left side variable appears also on the right hand side)

That's an oversimplification. There are other cases in addition to just
the left side variable appearing on the right hand side. The definition
of the standard is that the right-hand side is evaluated first. Anything
else is an optimzation. The optimization might make things faster,but it
does not change the meaning (not if it is a valid optimization).

And no, the compiler does not "convert matrix operations into do loops
and then compile." There have been cases of preprocessors that do things
like that, but that's not an accurate description of how compilers work
on these things. If Andy said that, he was trying to make a simple
analogy that he thougt you might understand.

And by the way, no the loop you showed is not a correct equivalent in
all cases. In particular, it does not catch the many-to-one assignment
issue that can come up with vector subscripts. The loop is valid even if
several of the itabsrt(1:8,2) elements have the same value; the array
assignment is not.

But you miss dpb's point anyway. It is a very important point, Missing
it will cause you make lots of errors eventually. You do *NOT* look at
context when figuring out what an expression means. In this case, you do
not look at the left-hand side to figure out what the right-hand side
means. You first figure out what the right-hand side means as it stands
alone. If you want to *LATER* think about how the compiler might
optimize things, it is really important that you consider that only
later. This is more than just some abstract conceptual thing that you
don't have to worry about. If you think about it any other way, you will
make errors - almost guaranteed.

It applies to more than just arrays. A quite common case is for people
to write something like

  some_double_precision_variable = 0.1

and think that this will give them a double precision aproximation to
0.1. It won't (except with some compilers that try to guess what you
meant instead of what you said, sometimes guessing incorrectly). The
right-hand side is a single precision approximation of 0.1. That's
because you don't look at the left-hand side to figure out what the
right-hand side means. The assignment does not gain you back the
precision that was lost in making the single-precision approximation on
the right-hand side.

-- 
Richard Maine                    | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle           |  -- Mark Twain

On Oct 5, 5:08=A0pm, nos...@see.signature (Richard Maine) wrote:
> The Star King <j...@npl.co.uk> wrote:
>
> > > The rhs is evaluated _COMPLETELY_in_ignorance_ of the LHS and that
> > > result is then _ASSIGNED_ to the LHS.
>
> > This may be conceptually how one should think of it. But what the
> > compiler may actually do is convert matrix operations into do loops
> > and then compile. This is just how g95 works, I asked Andy Vaught
> > about it. That is, when confronted with
>
> > ul(1:8,:)=3Duu(3,itabsrt(1:8,2),4,:)
>
> > g95 will not create temporary arrays. (it will create temporaries if
> > the left side variable appears also on the right hand side)
>
> That's an oversimplification. There are other cases in addition to just
> the left side variable appearing on the right hand side. The definition
> of the standard is that the right-hand side is evaluated first. Anything
> else is an optimzation. The optimization might make things faster,but it
> does not change the meaning (not if it is a valid optimization).
>
> And no, the compiler does not "convert matrix operations into do loops
> and then compile." There have been cases of preprocessors that do things
> like that, but that's not an accurate description of how compilers work
> on these things. If Andy said that, he was trying to make a simple
> analogy that he thougt you might understand.
>
> And by the way, no the loop you showed is not a correct equivalent in
> all cases. In particular, it does not catch the many-to-one assignment
> issue that can come up with vector subscripts. The loop is valid even if
> several of the itabsrt(1:8,2) elements have the same value; the array
> assignment is not.
>
> But you miss dpb's point anyway. It is a very important point, Missing
> it will cause you make lots of errors eventually. You do *NOT* look at
> context when figuring out what an expression means. In this case, you do
> not look at the left-hand side to figure out what the right-hand side
> means. You first figure out what the right-hand side means as it stands
> alone. If you want to *LATER* think about how the compiler might
> optimize things, it is really important that you consider that only
> later. This is more than just some abstract conceptual thing that you
> don't have to worry about. If you think about it any other way, you will
> make errors - almost guaranteed.
>
> It applies to more than just arrays. A quite common case is for people
> to write something like
>
> =A0 some_double_precision_variable =3D 0.1
>
> and think that this will give them a double precision aproximation to
> 0.1. It won't (except with some compilers that try to guess what you
> meant instead of what you said, sometimes guessing incorrectly). The
> right-hand side is a single precision approximation of 0.1. That's
> because you don't look at the left-hand side to figure out what the
> right-hand side means. The assignment does not gain you back the
> precision that was lost in making the single-precision approximation on
> the right-hand side.
>
> --
> Richard Maine =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0| Good judgment come=
s from experience;
> email: last name at domain . net | experience comes from bad judgment.
> domain: summertriangle =A0 =A0 =A0 =A0 =A0 | =A0-- Mark Twain

Thanks for your answer. I take your point about the many-to-one issue.
I guess Andy wouldn't mind if I copy what he said (after all g95 is
open source!)

Q:
> > Internally, does g95 convert matrix operations into do-loops before
> > compiling? It seems to me this would be quite a simple way
> of writing
> > a Fortran 90 compiler based on an existing F77 compiler.

A:
> There is a scalarization pass inside of the compiler that
> converts vector
> expressions into loops over scalar expressions.  Does that answer your
> questions?

The reason I asked him was I wanted to know whether it is "better" (ie
faster code) to write do-loops or array operations (personally, i
prefer do-loops, I can see what is going on more easily). Andy said
that, when equivalent, the two produce identical machine code.

The Star King wrote:
....
> Thanks for your answer. I take your point about the many-to-one issue.
> I guess Andy wouldn't mind if I copy what he said (after all g95 is
> open source!)
> 
> Q:
>>> Internally, does g95 convert matrix operations into do-loops before
>>> compiling? It seems to me this would be quite a simple way
>> of writing
>>> a Fortran 90 compiler based on an existing F77 compiler.
> 
> A:
>> There is a scalarization pass inside of the compiler that
>> converts vector
>> expressions into loops over scalar expressions.  Does that answer your
>> questions?
> 
> The reason I asked him was I wanted to know whether it is "better" (ie
> faster code) to write do-loops or array operations (personally, i
> prefer do-loops, I can see what is going on more easily). Andy said
> that, when equivalent, the two produce identical machine code.

IMO you're over-generalizing the answer that Andy gave and 
underestimating the importance of Richard's amplification of my response 
as well in consideration of only a single subset of the general problem 
for a specific compiler... :(

--

The Star King <jfb@npl.co.uk> wrote:

> Thanks for your answer. I take your point about the many-to-one issue.
> I guess Andy wouldn't mind if I copy what he said (after all g95 is
> open source!)
> 
> Q:
> > > Internally, does g95 convert matrix operations into do-loops before
> > > compiling? It seems to me this would be quite a simple way
> > of writing
> > > a Fortran 90 compiler based on an existing F77 compiler.
> 
> A:
> > There is a scalarization pass inside of the compiler that
> > converts vector
> > expressions into loops over scalar expressions.  Does that answer your
> > questions?
> 
> The reason I asked him was I wanted to know whether it is "better" (ie
> faster code) to write do-loops or array operations (personally, i
> prefer do-loops, I can see what is going on more easily). Andy said
> that, when equivalent, the two produce identical machine code.

Yes, Andy's answer sounds much more sensible (not to surprising, as he
would certainly know). Note that he talks about this happening during a
"scalarization pass inside of the compiler". That is a lot different
from your "before compiling". Similarly, Andy talks about loops, but not
"DO loops" as you did. A loop is a general concept. A DO loop is a piece
Fortran source code. When Andy says that scalarization produces loops,
that does not mean that it preprocesses the source code into Fortran DO
loops and then compiles that.

I also note an implication in your question to Andy. G95 was not "based
on an existing F77 compiler." In the relatively early days of f90, there
was a compiler from.... Parasoft I think it was (long enough ago that I
might have gotten it mixed up; I know it was one of those whose name
began with P; not PGI; might have been PSR, but I don't think so)....
that did work as you describe by preprocessing f90 source code into f77
and then compiling it with an existing f77 compiler. That compiler was a
flop. I bought a copy once (at work - not personally). I found it
completely unuseable. It was very buggy (even on f77-style code). And
even if it had been bug free, the implementation approach for modules
would have made it impractical to use for anything but toy applications.

-- 
Richard Maine                    | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle           |  -- Mark Twain

The Star King <jfb@npl.co.uk> wrote:
< On Oct 5, 5:08?pm, nos...@see.signature (Richard Maine) wrote:
(snip)

<> And no, the compiler does not "convert matrix operations into do loops
<> and then compile." There have been cases of preprocessors that do things
<> like that, but that's not an accurate description of how compilers work
<> on these things. If Andy said that, he was trying to make a simple
<> analogy that he thougt you might understand.

(snip)

< Thanks for your answer. I take your point about the many-to-one issue.
< I guess Andy wouldn't mind if I copy what he said (after all g95 is
< open source!)
 
< Q:
<> > Internally, does g95 convert matrix operations into do-loops 
<> > before compiling? It seems to me this would be quite a simple
<> > way of writing a Fortran 90 compiler based on an existing F77 
<> > compiler.
 
< A:
<> There is a scalarization pass inside of the compiler that
<> converts vector expressions into loops over scalar expressions.  
<> Does that answer your questions?

Most compilers are 'multiple pass', such that they go through
the source (or its internal representation) more than once.

In the early (and also not so early) days each pass wrote out an
intermediate file to be read by the next pass.  The combination
of all the passes is the compiler.  There are some cases, though,
where a compiler is used for some of the passes.  Before there
was g77, the popular unix Fortran compiler consisted of f2c
followed by a C compiler.  It is then not so obvious if one
should call f2c the first pass?  If a Fortran 90 version of
f2c existed, it would even expand array expressions as you say.
(C doesn't have array expressions like that.)
 
< The reason I asked him was I wanted to know whether it is "better" (ie
< faster code) to write do-loops or array operations (personally, i
< prefer do-loops, I can see what is going on more easily). Andy said
< that, when equivalent, the two produce identical machine code.

I believe my past advice on this is still true, though improvements
in optimization are ongoing.  It seems to me that for simpler
array operations the array expression is at least as fast.  For more
complicated ones it can be slower.  

In the case of a simple array assignment:

a(:,:,:,:,:)=b(:,:,:,:,:)

(with or without the colons, which here indicate the rank)

the compiler may be able to do it with one loop, assuming that it
knows the arrays are in contiguous storage.  On some processors,
the actual move could be done in one instruction!  

With more complicated expressions, temporary arrays may be used,
even when you know that they aren't needed.  There are many tricks
that have been done with the array intrinsic functions, many of
which do much more work than simple DO loops would have done.

-- glen

Richard Maine <nospam@see.signature> wrote:
< The Star King <jfb@npl.co.uk> wrote:
(snip)
 
<> The reason I asked him was I wanted to know whether it is "better" (ie
<> faster code) to write do-loops or array operations (personally, i
<> prefer do-loops, I can see what is going on more easily). Andy said
<> that, when equivalent, the two produce identical machine code.
 
< Yes, Andy's answer sounds much more sensible (not to surprising, as he
< would certainly know). Note that he talks about this happening during a
< "scalarization pass inside of the compiler". That is a lot different
< from your "before compiling". Similarly, Andy talks about loops, but not
< "DO loops" as you did. A loop is a general concept. A DO loop is a piece
< Fortran source code. When Andy says that scalarization produces loops,
< that does not mean that it preprocesses the source code into Fortran DO
< loops and then compiles that.

It might be, though, that the intermediate code generated is the
same as the intermediate code generated by DO loops.  It might even
be that the symbolic representation uses the symbol DO.
 
< I also note an implication in your question to Andy. G95 was not "based
< on an existing F77 compiler." In the relatively early days of f90, there
< was a compiler from.... 
(snip)
< that did work as you describe by preprocessing f90 source code into f77
< and then compiling it with an existing f77 compiler. 

Uncommon but possible, as you say.  The early C++ compilers were
preprecessors generating C code, followed by a C compiler.

Not so unusual, though, is to use the 'back end' (compiler terminology
for the later passes) of an existing compiler and to write a new
'front end' (the earlier passes).  

A common way of writing compilers first converts the input into
an intermediate code that is somewhat machine independent, but
with operations that most machines know how to do.  Operations
like MOVE and ADD.  The later, machine dependent, passes convert
that to specific instructions for the target machine.  When porting
to a new machine, only the later passes need to be rewritten.

Someone can say if I am wrong, but as I understand it the back
end for gfortran and g95 are pretty much the back end for the
gcc C compiler.  Among other results of doing that is the strange
code for routines with ENTRY statements.  (Visible in the debugger.)

Any routine using ENTRY compiles as one function (C terminology)
with all the possible dummy arguments (Fortran terminology), 
along with little functions for each entry point that call the
big one with the appropriate arguments.  (Most likely with
NULL pointers for the others.)  Since C doesn't have ENTRY,
the intermediate code doesn't have a way to describe it.

-- glen

On Oct 5, 2:49=A0pm, The Star King <j...@npl.co.uk> wrote:
> > > The intention is
>
> > > do i=3D1,8
> > > do j=3D1,3
> > > =A0 =A0 ul(i,j)=3Duu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j)
> > > enddo
> > > enddo
>
> > > It's probably not allowed but I can't think of a specific reason.
>
> > Now I'm confused with the rank of 'uu'... still, assuming 'uu' is a
> > rank-4 and that you forgot to declare one of the dimensions, all of
> > what I said above is still valid. You cannot associate a rank-2 array
> > on the left with a rank-4 array on the right! So your nested DO loop
> > is also not legal!
>
>
> I'm sure my do loops, at least, are legal!
>
> do i=3D1,8
> do j=3D1,3
> =A0 =A0 ul(i,j)=3Duu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j)
> enddo
> enddo
>
> itabsrt(i,1) is just a scalar integer. I can use it to index uu.
>
> uu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),j) is a scalar real. I am
> setting ul(i,j) (a scalar) to this number

Yes, you're right, I'm sorry! I was skipping lunch and not thinking
clearly then. ;)
Yes, the loop is legal, what you have inside is just a scalar-scalar
association. But you'll always need the DO loops over 'i' to make the
association; not on 'j', though. This is equivalent to your loop:

do i=3D1,8
   ul(i,:)=3Duu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),:)
enddo

And as long as your example goes, this is how far you can get with
colon matching. Since the 'i' is repeated in different dimensions, you
cannot translate that into colons. One simple example is the
calculation of the trace of a matrix (rank-2 array):

real :: A(10,10), trace

trace =3D 0.
do i=3D1,10
   trace =3D trace + A(i,i)
end do

To calculate the sum of diagonal components, you cannot get rid of the
loop and replace it with colons. SUM(A(:,:)), for example, will sum
all components of A; just because A has two identical colons in
different dimensions, it doesn't mean that the program will sum
components along matching indices. Same goes for your DO loop: you'll
always need the 'i' loop, in your particular case.

Best,

-- helvio

helvio <helvio.vairinhos@googlemail.com> wrote:
(snip)

< Yes, the loop is legal, what you have inside is just a scalar-scalar
< association. But you'll always need the DO loops over 'i' to make the
< association; not on 'j', though. This is equivalent to your loop:
 
< do i=1,8
<   ul(i,:)=uu(itabsrt(i,1),itabsrt(i,2),itabsrt(i,3),:)
< enddo
 
< And as long as your example goes, this is how far you can get with
< colon matching. Since the 'i' is repeated in different dimensions, you
< cannot translate that into colons. One simple example is the
< calculation of the trace of a matrix (rank-2 array):
 
< real :: A(10,10), trace
 
< trace = 0.
< do i=1,10
<   trace = trace + A(i,i)
< end do
 
< To calculate the sum of diagonal components, you cannot get rid of the
< loop and replace it with colons. SUM(A(:,:)), for example, will sum
< all components of A; just because A has two identical colons in
< different dimensions, it doesn't mean that the program will sum
< components along matching indices. Same goes for your DO loop: you'll
< always need the 'i' loop, in your particular case.

With the appropriate MASK argument to SUM it should work, though 
likely slower than the above DO loop.  Some people like to see what
can be done with one line of code, no matter how slow or unreadable
the result.

      trace=sum(a,reshape((/(((i.eq.j),i=1,10),j=1,10)/),shape(a)))

-- glen

On 2009-10-05, glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> Richard Maine <nospam@see.signature> wrote:
>< The Star King <jfb@npl.co.uk> wrote:
> (snip)
>  
><> The reason I asked him was I wanted to know whether it is "better" (ie
><> faster code) to write do-loops or array operations (personally, i
><> prefer do-loops, I can see what is going on more easily). Andy said
><> that, when equivalent, the two produce identical machine code.
>  
>< Yes, Andy's answer sounds much more sensible (not to surprising, as he
>< would certainly know). Note that he talks about this happening during a
>< "scalarization pass inside of the compiler". That is a lot different
>< from your "before compiling". Similarly, Andy talks about loops, but not
>< "DO loops" as you did. A loop is a general concept. A DO loop is a piece
>< Fortran source code. When Andy says that scalarization produces loops,
>< that does not mean that it preprocesses the source code into Fortran DO
>< loops and then compiles that.
>
> It might be, though, that the intermediate code generated is the
> same as the intermediate code generated by DO loops.  It might even
> be that the symbolic representation uses the symbol DO.

No, it's not called "DO"; you can check for yourself with
-fdump-tree-original.

> Not so unusual, though, is to use the 'back end' (compiler terminology
> for the later passes) of an existing compiler and to write a new
> 'front end' (the earlier passes).

Nowadays the "standard" organization for a production quality compiler
tends to have somewhat separated frontend (parsing, lexing, converting
to some more or less language-agnostic tree description), middle-end
(optimization passes working on trees while progressively "lowering"
the trees), and finally the back-end (create target specific
assembler along with some target-specific optimizations).

> Someone can say if I am wrong, but as I understand it the back
> end for gfortran and g95 are pretty much the back end for the
> gcc C compiler.

Not only the back end, but also the middle end is shared among all the
frontends.

> Among other results of doing that is the strange
> code for routines with ENTRY statements.  (Visible in the debugger.)
>
> Any routine using ENTRY compiles as one function (C terminology)
> with all the possible dummy arguments (Fortran terminology), 
> along with little functions for each entry point that call the
> big one with the appropriate arguments.  (Most likely with
> NULL pointers for the others.)  Since C doesn't have ENTRY,
> the intermediate code doesn't have a way to describe it.

Well yes, the intermediate representation that the front-ends hand
over to the middle end was designed by C programmers, and looks
somewhat like "pseudo-C". Specifically wrt ENTRY, another problem IIRC
is that unix ABI's and linkers don't tend to support jumping into the
middle of a function, and adding that would presumably be quite a lot
of effort, for little gain. However, the debugging issue you mention
could presumably be fixed by fixing the debug info generation, while
keeping the actual implementation as-is.


-- 
JB