[Nek5000-users] reshape array

[nek5000-users at lists.mcs.anl.gov]

PS - please note that the addressing formula should read:


[ z(i,j,k,e) ] = [z(1,1,1,1)] +
                          (i-1)+(j-1)*lx1+(k-1)*lx1*ly1+(e-1)*lx1*ly1*lz1)*wdsize

________________________________
From: Nek5000-users <nek5000-users-bounces at lists.mcs.anl.gov> on behalf of nek5000-users at lists.mcs.anl.gov <nek5000-users at lists.mcs.anl.gov>
Sent: Wednesday, August 15, 2018 9:34:29 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: Re: [Nek5000-users] reshape array



f77 array declarations do either one or two things, depending on

whether the array is passed as an argument or is declared within

the scope of the routine or function in question.


If the array is in a common block or is declared for the first time

within a subroutine/function, then the statement sets aside memory.


For example:


        subroutine blah(z)


        include 'SIZE'

        real x(lx1,ly1,lz1,lelt)

        common /c1/ y(lx1,ly1,lz1,lelt)


        real z(lx1,ly1,lz1,lelt)


        integer e




results in two arrays, x and y, be declared to hold n words of memory,

where n = lx1*ly1*lz1*lelt


The line "real z(...)" does not declare any memory.


All three statements indicate how to reference an entry in any

of the arrays by the formula:



            [ z(i,j,k,e) ] = [z(1,1,1,1)] +

                                  (i+(j-1)*lx1+(k-1)*lx1*ly1+(e-1)*lx1*ly1*lz1)*wdsize


where [.] implies the address of the argument.   This formula applies

equally to x,y, and z in the above example.  f77 does not care about

i,j,k, or e (save that they must be integer).  It assumes you know what

you are doing when you compute an address.   Thus the following are

almost equivalent:



      do e=1,nelt

      do k=1,lz1

      do j=1,ly1

      do i=1,lx1

           x(i,j,k,e)=z(i,j,k,e)

      enddo

      enddo

      enddo

      enddo


      n = lx1*ly1*lz1*nelt   (Note: nelt, not lelt)


      do i=1,n

           x(i,1,1,1) = z(i,1,1,1)

      enddo


The latter, however, is _much_ better than the former.  Why?


Most compilers will not vectorize more than one loop deep (the "i" loop

in the first of the two preceding examples).   You thus have a relatively short loop for vectorization; whereas the second, single-loop, form will vectorize.  Moreover, the second form makes the code much easier to use.   The essential point is that x and z are contiguous arrays --- which is the case throughout nek.


NOTE:  many users are hesitant about using the 2nd form, despite the fact that it is vastly superior to the first.  There is no need to be uneasy about this form, however.  It has been used in Nek for over 30 years on hundreds of the world's fastest computers.  It is a well-know fortran trick for high performance.   Please look through the Nek5000 source for more examples of this type.



________________________________
From: Nek5000-users <nek5000-users-bounces at lists.mcs.anl.gov> on behalf of nek5000-users at lists.mcs.anl.gov <nek5000-users at lists.mcs.anl.gov>
Sent: Wednesday, August 15, 2018 4:16:15 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: [Nek5000-users] reshape array


Hi all,



All average streamwise velocity are stored in this array



“common /avg/ uavg(lx1,ly1,lz1,lelv)”



but I want this array to be reshaped into



“common /abcd/ abcd(lx1*ly1*lz1,lelv)”



The following didn’t work:





common /avg/   uavg(lx1,ly1,lz1,lelv)

common /abcd/  abcd(lx1*ly1*lz1,lelv)





abcd=reshape(uavg,(/ lx1*ly1*lz1,lelv /))





This is how array reshaping works in fortran. Unfortunately, didn’t work in nek. Anyone knows why?



Many thanks in advance for your replies.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180815/a0ef3c4b/attachment.html>