Fortran - Allocatable derived type with allocatable components in fortran 90

Author

Allocatable derived type with allocatable components in fortran 90

2006-03-22, 7:03 pm

I want to do something similar to what I show below. Needless to say,
what I have shown below will fail miserably. But can someone tell me how
to achieve what I want to achieve ?

integer :: N_dimensions
integer :: N_protons

type particle
real, dimension(:), pointer :: position
real, dimension(:), pointer :: velocity
end type particle

type (particle), dimension (:), allocatable :: proton

allocate (proton%position(N_dimensions))
allocate (proton%velocity(N_dimensions))

allocate (proton(N_protons))

Thanks in advance,
Shriram

Richard E Maine

2006-03-22, 7:03 pm

<dont@spam.me> wrote:

> I want to do something similar to what I show below. Needless to say,
> what I have shown below will fail miserably. But can someone tell me how
> to achieve what I want to achieve ?

Well, I notethat your subject line doesn't quite match. What you have
here is not a derived type with allocatable cmponents. It is a derived
type with pointer components. Pointer components can often be used to
"fake" allocatable ones, but it results in non-intuitive behavior in
some cases.

If your compilers support allocatable components (i.e. the allocatable
TR to f95), you might be advised to use that. To do so, just change
pointer to allocatable in the component declarations. All the rest will
be the same.

> integer :: N_dimensions
> integer :: N_protons
>
> type particle
> real, dimension(:), pointer :: position
> real, dimension(:), pointer :: velocity
> end type particle
>
> type (particle), dimension (:), allocatable :: proton
>
> allocate (proton%position(N_dimensions))
> allocate (proton%velocity(N_dimensions))
>
> allocate (proton(N_protons))

Your problem here is mostly that your allocations are backwards. You
need to allocate proton first, before doing anything with any of its
components. That includes allocating the components. Also, you'll need
to allocate the components of each proton separately. For example

allocate(proton(N_protons))
do i = 1 , N_protons
allocate(proton(i)%position(N_dimensions
))
allocate(protons(i)%velocity(N_dimension
s))
end do

P.S. Is N_dimensions anything other that a constant 3? If it is always
3, then things can be simpler. If, on the other hand, your protons might
be moving in 11-dimensional (or whatever the current favorite number is)
string space, I suppose you do need the allocatable bit for the
components.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain| experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Richard E Maine

2006-03-23, 7:04 pm

Shriram Ramanathan <ramanath@cems.umn.edu> wrote:

> Richard and Joe, Thanks for your answers. One more question
>
> Does the allocation method mentioned by Richard ensure that all data for
> a given particle (say the 112th proton) is stored contiguously in memory.
> I just want to ensure that when I access say the position of the
> 112th proton, then the velocity associated with the 112th proton is also
> present in the cache (or rather there is a high probability that it is in
> the cache). That way, cache misses can be greatly reduced.
>
> If this method does not ensure that particle data is stored contiguously,
> then what is the best way to handle this ?

Well, I mentioned a couple of methods. I spelled out one in code, but I
mentioned several. One of them does - that's the method (also mentioned
by Joe, with sample syntax shown) of making the component dimensions
fixed size (presumably 3). I didn't hear an answer to the question of
whether the fixed size would work for you. If a fixed size does work,
then I'd strongly recommend it. It is fundamentally simpler. And yes, it
guarantees contiguity (not quite in so many words in the standard, but
close, and certainly in practice).

The methods that involve allocatable or pointer components do not
guarantee that. Indeed, they pretty much guarantee the opposite. The
actual data array for an allocatable or pointer component is going to be
at some other location in memory, separate from that of any other
components. One could manage to force contiguity with pointers, but it
would be a mess (allocate a large target array elswehere and make all
the pointers point into areas of it that you parcel out); you'd have to
have really good reason to do that. And you'd have good chance of
hurting rather than helping performance because of pointer aliasing
isssues.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain| experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

s8ngsu3@yahoo.com

2006-03-24, 8:01 am

ramanath@cems.umn.edu wrote:

> Just in case someone has ideas on how to do this one, let me
> mention what I want to do. Briefly, I will be having a simulation domain
> with N_particles. This simulation domain will be divided into
> N_cells with each cell containing (N_particles_cell) particles. Most
> of the steps in the simulation would be limited to operating on particles
> residing in a cell. # of Operations involving inter-cell communication
> between data in two different cells is much smaller. So, I am trying to
> fit all data related to a cell within cache.
>
> I basically wanted to define a derived type like "cell_data" which would
> contain all data related to the cell
>
> type cell_data
> integer :: cell_volume
> integer :: No_of_particles_in_cell ! (which varires as the simulation
> proceeds)
> integer, dimension(:), allocatable :: list_of_particles_in_cell
> ! size of array will be (max_allowed_no_of_particles_in_cell)
> real, dimension(:,:), allocatable :: position_particles
> ! size of array will be (no_of_dimensions, max_allowed_no_of_particles_in_cell)
>
> real, dimension(:,:) :: velocity_particles ! same size as above
> end type cell_data
>
> type(cell_data), dimension(:) :: cells
> allocate (cells(N_cells))

I guess you want to create the so-called Verlet list, in order to
identify quickly the particles that interact at the given time moment.
In this case, you may want to store only the particle numbers in the
cells, and keep the other data (coordinates, velocities, accelerations,
etc.) in simple arrays. Having the particle numbers in the cells, you
can access this data. The advantage is that if you solve the equations
of motion, you can use the whole-array operations, which is fast; e.g,
if you want to update the positions:

real, allocatable:: positions(:,:), velocities(:,:), ...
....
positions=positions + dt * velocities

smth like that.

As far as I understand, you expect to speed up by ensuring that the
relevant portions of these arrays fit into the CPU cache? But usually,
when solving the equations of motion, you have to update the data of
all the particles during a time step, therefore, you have to operate on
the entire arrays anyway.

Or did I misunderstood something?

Sponsored Links