Aleksandar Donev wrote:

> In the lack of F95 (or at present F2x) features to have such references,
> I have been trying today to emulate this, and am stuck with a question
> that I hope those experienced in array sections might be able to answer.
>
> Let's start with an actual example--weekly weather records over a 2D
> cartesian spatial grid of stations. We have the following declarations:
>
>          TYPE :: Hourly_Record
>             REAL (KIND=r_wp) :: temperature (3) = 0.0
>             LOGICAL (KIND=l_byte) :: synny = .TRUE.
>          END TYPE
>
>          TYPE :: Daily_Record
>             TYPE (Hourly_Record), DIMENSION (24) :: hourly_records
>             INTEGER (KIND=i_sp) :: sunrise = 7, sunset = 18
>          END TYPE
>
>          TYPE :: Weekly_Record
>             TYPE (Daily_Record), DIMENSION (7) :: daily_records
>             REAL (KIND=r_sp) :: forecast_success (5)
>          END TYPE
>
>          INTEGER, PARAMETER :: n_x = 10, n_y = 15
>          TYPE (Weekly_Record), DIMENSION (n_x, n_y), TARGET ::
> weekly_records ! The grid of records
>
> The point of the discussion is that one would like to associate a rank-5
> array pointer with the array section of all temperatures, but at present
> one cannot do this. So what I want to achieve is:
>
> REAL(KIND=r_wp), DIMENSION(:,:,:,:,:), POINTER :: temperatures
>
> temperatures=>weekly_records[:,:]%daily_records[:]%hourly_records[:]%temperature[:]
> ! Not legal at present

You are asking to have a data structure supporting dual views to users, as
a multi-level structure and as a flat array.  The multi-level structure view
is defined as the storage, while the flat array view is used as an alias.

IMHO, a realistic way in Fortran should really be in the opposite.  The flat
array should be the primary view as the storage, while the multi-level
structure view should be the alias.

This is what I mean:

> TYPE :: Hourly_Record
     REAL (KIND=r_wp),pointer,dimension(:) :: temperature ! <- this is different
>    LOGICAL (KIND=l_byte) :: synny = .TRUE.
> END TYPE
>
> TYPE :: Daily_Record
>    TYPE (Hourly_Record), DIMENSION (24) :: hourly_records
>    INTEGER (KIND=i_sp) :: sunrise = 7, sunset = 18
> END TYPE
>
> TYPE :: Weekly_Record
>    TYPE (Daily_Record), DIMENSION (7) :: daily_records
>    REAL (KIND=r_sp) :: forecast_success (5)
> END TYPE

  type WorldWeeklyRecords
     REAL (KIND=r_wp),pointer,dimension(:,:,:,:,:) :: storage
     type(Weekly_Record),pointer,dimension(:,:) :: weekly_records
  end WorldWeeklyRecords

In an application, the storage is allocated first.

> INTEGER, PARAMETER :: n_x = 10, n_y = 15
  type(WorldWeekly_Record) :: worldweekly_records

  allocate(worldweekly_records%storage(3,24,7,n_x,n_y))
  allocate(worldweekly_records%weekly_records(n_x,n_y))

Then the loops aliasing the components of weekly to alltemp%storage must
be done:

  do j=1,n_y; do i=1,n_x; do id=1,7; do ih=1,24
    worldweekly_records%weekly(i,j)%daily(id)%hourly(ih)%temperature => &
    worldweekly_records%storage(:,ih,id,i,j)
  end do; end do; end do

Then you fly.

This approach is obviously not straight forward, but should work nicely.

>
> which should produce an array pointer of shape [3,24,7,10,15]. Then one
> can take subsections of this pointer to for example find the
> temperatures on odd days of the week in the morning hours along a given
> meridian.

By creating necessary module procedures of pointer types and maybe additional
components (pointers) for the components of WorldWeekly_Record, one can
almost dice the data in anyway applications need.

If one does it right, the code could also be highly reusable.

>
> I have been trying to make some routines that will create this pointer
> for me, but I am having trouble.
>
> Assume I know the following (I find this using C functions for example)
> about the actual memory layout of the weekly_records structure:
>
> The type Hourly_Records take 4 words. I'll denote HR=4*1
> There are 24 of these in a daily record + 2 trailing words=96+2=98
> words. DR=24*HR+2
> Then there are 7 of these in a weekly record + 5 words=686+5=691 words.
> WR=7*DR+5
> Finally, there are 10 of these per column of the grid=69100 words.
> COL=10*WR
> And there are 15 of these total, i.e. 1036500 words. TOTAL=15*COL
>
> So the whole memory layout of the datastructure is
> 15*(10*(7*(24*(4)+2)+5)) in words.
>
> Now, I claimed, and others seem to agree, that the section
> weekly_records[:,:]%daily_records[:]%hourly_records[:]%temperature[:]
> can indeed be described as a *regular* rank-5 array section. HOW?!?

That is all nice, but that is the old way.  I am aware that many people
are thinking this way these days.  But neo-Fortran programmer should not
(have to) do that anymore.  :-)

Likely, this won't solve your problem, since you might not have any
freedom to redefine the source data already in a multi-leveled structure.
Even so, I think it is still important to rise the awareness of Fortran
friendly designs.

The old FORTRAN prefers arrays.  The new Fortran still prefers arrays.
Any ADT is an overhead on top of arrays.  The other way arround just
won't be natural, atleast for the time being.

Jing

>
> I would say (but may be wrong!) that *any* array section of rank r can
> be described in two steps:
>
> 1) First give an underlying (call it base) *contiguous* array of rank r
> (the memory buffer of 1036500 words, shaped as a rank-r array)
> 2) Then take an appropriate Fortran section slice through this buffer,
> base(lb1:ub1:stride1,...).
>
> The problem I am facing is finding the shape of the base array. If I
> could do this, then for example I could make an array pointer to the
> base buffer (using my emulation of C_F_POINTER from F2x) and use a
> section of it to make the desired pointer temperatures.
>
> But in the example above I cannot see how to describe the desired
> section in the form base(section). Does someone see how to do this? Or
> am I thinking the wrong way here?
>
> Any help with this indexing issue will be appreciated. I will of course
> post the final code for public use.
>
> Thanks,
> Aleksandar
>
> --
> __________________________________
> Aleksandar Donev
> Complex Materials Theory Group (http://cherrypit.princeton.edu/)
> Princeton Materials Institute & Program in Applied and Computational Mathematics
> @ Princeton University
> Address:
>    419 Bowen Hall, 70 Prospect Avenue
>    Princeton University
>    Princeton, NJ 08540-5211
> E-mail: [log in to unmask]
> WWW: http://atom.princeton.edu/donev
> Phone: (609) 258-2775
> Fax: (609) 258-6878
> __________________________________

--
________________________________ _-__-_-_ _-___--- _____________________
Jing Guo, [log in to unmask], (301)614-6172(o), (301)614-6297(fx)
Data Assimilation Office, Code 910.3, NASA/GSFC, Greenbelt, MD 20771