How implement a pointer to allocatable array inside a derived data type?

Help
1

I have a mesh like so:

|  .  |   .   |   .   |  ...   |
0     1       2       3        nc

To store the position of the cell boundaries, centers, etc., I have created the following derived data type:

testgrid1.f90 (2.1 KB)

   type :: grid
        !! 1D grid
      character(20) :: name = "x [-]"
        !! variable name
      real(rk), allocatable :: edges(:)
        !! vector(0:nc) of cell edges
      real(rk), allocatable :: center(:)
        !! vector(nc) of cell centers, \( x_i \)
      real(rk), allocatable :: width(:)
        !! vector(nc) of cell widths,  \( x_{i+1/2} - x_{i-1/2} \)
      real(rk), allocatable :: left(:)
        !! vector(nc) of left cell boundaries, \( x_{i-1/2} \)
      real(rk), allocatable :: right(:)
        !! vector(nc) of right cell boundaries, , \( x_{i+1/2} \)
      integer :: ncells
        !! number of cells
   contains
      procedure, pass(self) :: new
   end type grid

This implementation works well and is convenient to use, but it is far from efficient in terms of memory, because the left=edges(0:nc-1) and right=edges(1:nc). The same info is being stored 3 times.
To “solve” this problem, I naively thought about redefining left and right to be pointers that would point to edges, i.e, something like:

   type :: grid
        !! 1D grid - DOES NOT WORK
      character(20) :: name = "x [-]"
        !! variable name
      real(rk), allocatable, TARGET:: edges(:)
        !! vector(0:nc) of cell edges
      real(rk), allocatable :: center(:)
        !! vector(nc) of cell centers, \( x_i \)
      real(rk), allocatable :: width(:)
        !! vector(nc) of cell widths,  \( x_{i+1/2} - x_{i-1/2} \)
      real(rk), dimension(:), POINTER :: left
        !! vector(nc) of left cell boundaries, \( x_{i-1/2} \)
      real(rk), dimension(:), POINTER :: right
        !! vector(nc) of right cell boundaries, , \( x_{i+1/2} \)
      integer :: ncells
        !! number of cells
   contains
      procedure, pass(self) :: new
   end type grid

But his kind of silly trick is completely illegal! :shushing_face: The attribute target is not allowed inside a derived data type.
So, what is proper way to achieve the intended result?
Thanks!

2

You don’t need to declare a target inside your derived type, rather the instance of the derived type needs to have the target attribute

subroutine new(self)
  class(grid), intent(out), target :: self
  ! ...
  self%left => self%edges(0:self%ncells - 1)
  self%right => self%edges(1:self%ncells)
end subroutine new
3

Thanks once again, @awvwgk !
It seems extremely obvious, now that you explained it to me!

4

FWIW, I think you could just define an elemental type-bound function in your module for them:

elemental real(rk) function left(self,celli)
   class(grid), intent(in) :: self
   integer, intent(in) :: celli
   left = self%edges(celli-1)
end function left

This will be inlined by any compilers, and lets you avoid the pointer thing at all.

5

Is this really true? How can it (practically) be inlined when the procedure is only dispatched at runtime due to being type-bound?

Tried it at Compiler Explorer for comparison of the pointer and elemental implementation.

6

Attention @awvwgk , @HugoMVale ,

For a standard-conforming program, note it is not adequate for the parameter of the type-bound procedure to have the TARGET attribute, the passed object too must have the TARGET attribute in order for the associations to persist beyond the end of the procedure invocation.

For example, in your reproducer. the gx object in the testgrid main program shall be declared with a TARGET attribute.

Otherwise, the pointer associations such as self%left => self%edges(.. only apply during the execution phase of the new procedure.

7

While technically true, I think the pointer will remain valid even outside the new procedure as the actual target is unlikely to move after leaving the procedure. To move the allocatable outside of new and keep the pointer valid the target attribute is required.

But sure, it can become an actual issue, therefore I try to avoid pointer components in derived types.

8

I’m not great at reading assembly, but in your example, it’s enough to turn on some non-zero level of optimization to achieve it: -O1 with ifort and -O2 with gfortran. Compiler experts will have better comments than mine, but I believe polymorphism only kicks in when you have a polymorphic variable (declared class(whatever), allocatable), while if your variable is statically typed (type(whatever_1)) the compiler knows exactly what functions have to be called, and can inline them

9

This may be true within the same file, but I don’t think you can make this statement more generally if this function could be called from a separately-compiled source file.
This is important because it could be the difference between a loop being vectorised or not.

Unless the compiler is storing implementation code in .mod files[1] (which I don’t think gfortran does) or you have specified link-time optimisation (lto), then the compiler has no way to inline the definition of a function or subroutine into a separately compiled file. This is one of the advantages of C++ in that included header files allow the compiler to see and inline the body of any functions that are declared and defined there.

By comparison, the simple pointer solution does not suffer from this limitation because accessing an element in the pointer array can be inlined and vectorised.

  1. The mod files produced by ifort are often larger than those produced by gfortran, so I wondered whether they store code definitions to allow inlining, but I haven’t checked this. ↩︎

10

Right, mileage varies with compiler when it comes to link-time optimization. That’s also fully supported by gcc/gfortran afaict.

Going back to the problem, of course the elemental solution is not going to be super-fast if that function is used alot from outside the module; my guess is most of the finite-difference work will be done from within the module itself.

Your exact comment is one of the reasons I support the addition of the inline keyword for Fortran btw:
https://github.com/j3-fortran/fortran_proposals/issues/229