Array of pointers in Fortran

Question

I have been writing a large Fortran program for thermodynamic calculations for almost 10 years and when I started I was new to the new Fortran standard (I was familiar with F77 and too old to learn something else). I found the new TYPE constructs very nice and have used them frequently but I was not aware of some limitations, such as it was not allowed to create arrays of pointers, I have discovered that later.

Now I am correcting some of my old code and I am surprised to find inside a record declaration for: TYPE gtp_phase_add

a declaration: TYPE(tpfun_expression), dimension(:), pointer :: explink

where explink is used to point to another structure containing a mathematical expression. This has not generated any compilation errors (I normally use gfortran but I have complied this program with intel fortran also). When I saw this old code (written some 10 year ago) I thought there is an "allocatable" missing but adding that gave a compilation error.

I made a minimal complete program to mimic how this is used:

MODULE test1
  implicit none
  
  TYPE tpfun_expression
     integer nc
     double precision, allocatable, dimension(:) :: coeffs
     integer, allocatable, dimension(:) :: powers
  END type tpfun_expression

  TYPE gtp_phase_add
!**************************************************************************
! My question is if it is correct Fortran to have an array of pointers here
     TYPE(tpfun_expression), dimension(:), pointer :: explink
!**************************************************************************
     TYPE(gtp_phase_add), pointer :: nextadd
  END TYPE gtp_phase_add

contains

  subroutine create_tpfun(n,coeffs,powers,exp)
    integer n,i,powers(*)
    double precision coeffs(*)
    type(tpfun_expression), pointer :: exp
    allocate(exp%coeffs(n))
    allocate(exp%powers(n))
    exp%nc=n
    do i=1,n
       exp%coeffs(i)=coeffs(i)
       exp%powers(i)=powers(i)
    enddo
    return
  end subroutine create_tpfun

  subroutine create_addrec(typ,this)
    integer typ,n,m
    TYPE(tpfun_expression), target :: exp1
    TYPE(tpfun_expression), pointer :: exp2
    TYPE(gtp_phase_add), pointer :: this
    integer ipow(4)
    double precision coeffs(4)
!
!**************************************************************************
! here I allocate a pointer array
    allocate(this%explink(typ))
!**************************************************************************
    if(typ.eq.1) then
       do m=1,4
          ipow(m)=m-1
          coeffs(m)=2.0D0*m
       enddo
       exp2=>this%explink(1)
       call create_tpfun(4,coeffs,ipow,exp2)
    else
       do m=1,4
          ipow(m)=m-1
          coeffs(m)=3.0D0
       enddo
       exp2=>this%explink(1)
       call create_tpfun(4,coeffs,ipow,exp2)
       do m=1,3
          ipow(m)=1-m
          coeffs(m)=5.0D0
       enddo
       exp2=>this%explink(2)
       call create_tpfun(3,coeffs,ipow,exp2)
    endif
    return
  end subroutine create_addrec

end MODULE test1

program main
  use test1
  integer n,m,j,k,q
  TYPE(gtp_phase_add), target :: addrec
  TYPE(gtp_phase_add), pointer :: next,first
  TYPE(tpfun_expression) :: exp
  
  first=>addrec
  next=>addrec
  write(*,*)'Creating addrec 1'
  call create_addrec(1,next)
  allocate(next%nextadd)
  write(*,*)'Creating addrec 2'
  next=>next%nextadd
  call create_addrec(2,next)
! just a check that the functions are correct
  write(*,*)'Listing functions in all addrecs'
  next=>first
  q=0
  do while(associated(next))
     q=q+1
     write(*,*)'Addition record ',q
     n=size(next%explink)
     do m=1,n
        k=next%explink(m)%nc
        write(*,10)(next%explink(m)%coeffs(j),next%explink(m)%powers(j),j=1,k)
     enddo
10   format(10(F6.3,1x,i3))
     next=>next%nextadd
  enddo
end program main

This works as I expect, I am only surprised that I allowed to declare an array of pointers so I would like to know if this is correct Fortran.' And sorry if I am not able to understand how to edit this more elegantly.

Answer 1

This is not an answer but a follow up to explain what was my confusion.

! Another test of fortran pointer arrays
  module example1
    implicit none
  
    type an_element
       character symbol*2
       type(an_element), pointer :: nextel
    end type an_element

    type subsystyp1
       character sysname*24
! here it seems I can have an array of pointers
       type(an_element), dimension(:), pointer :: ellista
    end type subsystyp1
    
    type elarray
! this is a type which contain just a pointer
       type(an_element), pointer :: p1
    end type elarray

    type subsystyp2
       character sysname*24
! here an array of a type which contain a pointer can be allocated
       type(elarray), dimension(:), allocatable :: ellista
    end type subsystyp2

  end module example1

  program test
    use example1
    integer ii
    type(an_element), target :: element
    type(an_element), pointer :: firstelem
    type(an_element), pointer :: nextelem
    type(subsystyp1) :: sys1
    type(subsystyp2) :: sys2
!
! create a linked list of elements
    element%symbol='Al'
    firstelem=>element
!
    allocate(element%nextel); nextelem=>element%nextel
    nextelem%symbol='Be'
    allocate(nextelem%nextel); nextelem=>nextelem%nextel
    nextelem%symbol='C'
    allocate(nextelem%nextel); nextelem=>nextelem%nextel
    nextelem%symbol='Fe'
    allocate(nextelem%nextel); nextelem=>nextelem%nextel
    nextelem%symbol='Mn'
    nullify(nextelem%nextel)
! list all the elements in the list
    nextelem=>firstelem
    write(*,10,advance='no')'All elements: '
    do while(associated(nextelem))
       write(*,10,advance='no')nextelem%symbol
       nextelem=>nextelem%nextel
10     format(a,1x)
    enddo
    write(*,*)
! Allocate a subsystem using subsystyp1  and set pointers to same elements
    sys1%sysname='ternary sys1'
    allocate(sys1%ellista(3))
! I cannot assign this as a pointer, that creates a compiler error
!    sys1%ellista(1)=>firstelem
    sys1%ellista(1)=firstelem
    sys1%ellista(2)=firstelem%nextel%nextel
    sys1%ellista(3)=sys1%ellista(2)%nextel
    write(*,20)trim(sys1%sysname),(sys1%ellista(ii)%symbol,ii=1,3)
20  format(a,': ',10(a,1x))
! Allocate a subsystem using subsystyp2 and set pointers to elements
    allocate(sys2%ellista(3))
    sys2%sysname='ternary sys2'
    sys2%ellista(1)%p1=>firstelem
    sys2%ellista(2)%p1=>firstelem%nextel%nextel
    sys2%ellista(3)%p1=>sys2%ellista(2)%p1%nextel
!
    write(*,20)trim(sys2%sysname),(sys2%ellista(ii)%p1%symbol,ii=1,3)
! So far so good ...
! Now change the element symbol in the one of the elements
    write(*,*)'Change the name of third element to "He" in main list'
    firstelem%nextel%nextel%symbol='He'
! list the whole list
    nextelem=>firstelem
    write(*,10,advance='no')'All elements: '
    do while(associated(nextelem))
       write(*,10,advance='no')nextelem%symbol
       nextelem=>nextelem%nextel
    enddo
    write(*,*)
! In sys1 the element name has nnot changed because ellista was assigned with "="
    write(*,20)trim(sys1%sysname),(sys1%ellista(ii)%symbol,ii=1,3)
! In sys2 the element name has changed because p1 is a pointer
    write(*,20)trim(sys2%sysname),(sys2%ellista(ii)%p1%symbol,ii=1,3)
!
  end program test

The output will be

All elements:  Al Be C  Fe Mn
ternary sys1: Al C  Fe
ternary sys2: Al C  Fe
 Change the name of third element to "He" in main list
All elements:  Al Be He Fe Mn
ternary sys1: Al C  Fe
ternary sys2: Al He Fe

A pointer assigned by "=" to a record will make a copy of that record. If later there are changes in the record that will not be included in the copy. But if a pointer is assigned by "=>" the changes will be available. That is trivial but if the assignment is not explicit, as by the subroutine call, the result can be confusing.

I know I am old fashioned and like to keep the "return" and I do not like the syntax of this implied do loop. To me it looks like an invention of a manager who pays his programmers by the number of lines of code they write.

Answer 2

@BoSundman: As you say in comments that you are still confused, I thought I would offer a frame challenge which might help clear some things up.

In particular, I recommend replacing your pointers to arrays with just using arrays. Your original use of pointers to arrays is completely valid Fortran, but I would argue that pointers to arrays are not the best tool for the job here.

I have re-written the code snippet you gave using what I hope is clearer Fortran:

module precision
  ! This is the modern way of defining floating point precision.
  integer, parameter :: dp = selected_real_kind(15,300)
end module

MODULE test1
  use precision
  implicit none
  
  type tpfun_expression
    ! 'nc' does not need to be stored separately, as nc=size(coeffs).
    real(dp), allocatable, dimension(:) :: coeffs
    integer,  allocatable, dimension(:) :: powers
  end type tpfun_expression

  type gtp_phase_add
    ! Rather than having 'explink' as a pointer to an array,
    !    consider using an allocatable array.
    type(tpfun_expression), allocatable, dimension(:) :: explink
    type(gtp_phase_add), pointer                      :: nextadd
  end type gtp_phase_add

contains

  ! Consider replacing subroutines with functions where appropriate.
  function create_tpfun(coeffs,powers) result(exp)
    ! Using 'intent(in)' tells the compiler (and anyone reading your code) that
    !    input arguments will not be modified.
    real(dp), intent(in)   :: coeffs(:)
    integer,  intent(in)   :: powers(:)
    ! The function returns a `type(tpfun_expression)` object, rather than
    !    working with a pointer.
    type(tpfun_expression) :: exp
    
    ! You can perform array copying as a single operation.
    exp%coeffs = coeffs
    exp%powers = powers
    
    ! There is no need for 'return' at the end of a function or subroutine.
    ! 'return' happens automatically when the end is reached.
  end function

  function create_addrec(typ) result(this)
    integer, intent(in) :: typ
    TYPE(gtp_phase_add) :: this
    
    integer,  allocatable :: ipow(:)
    real(dp), allocatable :: coeffs(:)
    
    integer :: m
    
    allocate(this%explink(typ))
    if(typ.eq.1) then
      ! The expression [(m-1,m=1,4)] uses an "implied do loop" to create an
      !    array in a single line of code. This is equivalent to the lines:
      !
      ! allocate(ipow(4))
      ! do m=1,4
      !   ipow(m) = m-1
      ! enddo
      !
      ipow = [(m-1,m=1,4)]
      coeffs = [(2.0_dp*m,m=1,4)]
      ! Since 'create_tpfun' is now a function,
      !    the syntax of calling it changes slightly.
      this%explink(1) = create_tpfun(coeffs,ipow)
    else
      ipow = [(m-1,m=1,4)]
      coeffs = [(3.0_dp,m=1,4)]
      this%explink(1) = create_tpfun(coeffs,ipow)
      
      ipow = [(1-m,m=1,3)]
      coeffs = [(5.0_dp,m=1,3)]
      this%explink(2) = create_tpfun(coeffs,ipow)
    endif
    
    ! The pointer 'this%nextadd' should be nullified.
    nullify(this%nextadd)
  end function
end module

program main
  use precision
  use test1
  implicit none
  
  TYPE(gtp_phase_add), target  :: first
  TYPE(gtp_phase_add), pointer :: next
  
  integer m,j,q
  
  write(*,*)'Creating addrec 1'
  next=>first
  next = create_addrec(1)
  write(*,*)'Creating addrec 2'
  allocate(next%nextadd)
  next=>next%nextadd
  next = create_addrec(2)
  
  ! just a check that the functions are correct
  write(*,*)'Listing functions in all addrecs'
  next=>first
  q=0
  do while(associated(next))
     q=q+1
     write(*,*)'Addition record ',q
     ! I have replaced 'n' and 'k' with 'size(next%explink)' and
     !    'size(next%explink(m)%coeffs)' respectively.
     do m=1,size(next%explink)
        ! I broke this line across multiple lines using '&',
        !    as it was a little long.
        write(*,10) ( next%explink(m)%coeffs(j),   &
                    & next%explink(m)%powers(j),   &
                    & j=1,                         &
                    & size(next%explink(m)%coeffs) )
     enddo
10   format(10(F6.3,1x,i3))
     next=>next%nextadd
  enddo
end program main

Answer 3

N.B. a variable integer, pointer :: ptr(:) is a pointer to an array (or an array slice), rather than an array of pointers. By this, I mean that this program is valid:

program test
  implicit none
  
  integer, target, allocatable :: foo(:)
  integer, pointer             :: ptr(:)
  
  foo = [1,2,3]
  
  ptr => foo
  write(*,*) ptr ! writes "1 2 3".
  
  ptr => foo(2:1:-1)
  write(*,*) ptr ! writes "2 1".
end program

but this program is not:

program test
  implicit none
  
  integer, target, allocatable :: foo(:)
  integer, target              :: bar
  integer, pointer             :: ptr(:)
  
  foo = [1,2,3]
  bar = 4
  
  ptr => foo
  ptr(2) => bar ! Not valid.
  write(*,*) ptr
end program

If you want something that acts as an array of pointers, you need to wrap the pointers in a class, for example:

program test
  implicit none
  
  type :: IntPointer
    integer, pointer :: ptr
  end type
  
  integer, target,  allocatable :: foo(:)
  integer, target               :: bar
  type(IntPointer), allocatable :: ptr(:)
  
  integer :: i
  
  foo = [1,2,3]
  bar = 4
  
  allocate(ptr(2))
  ptr(1)%ptr => foo(2)
  ptr(2)%ptr => bar
  
  write(*,*) (ptr(i)%ptr, i=1, size(ptr))
end program

This is also why a pointer cannot be allocatable: the pointer is not an array of separate pointers, each with their own target, but is rather a single pointer to an array (or array slice), plus some metadata about the size and shape of that array (or array slice). This means the memory requirements of the pointer are independent of the size() of the target, and so for memory management purposes an allocatable pointer wouldn't make sense.

When you call allocate(ptr) on a pointer :: ptr(:), you are not really allocating ptr, but rather allocating an (unnamed) array and then pointing ptr at that array.

Answer 4

Yes, you can have an array with the pointer attribute. Consider the trivial code

program foo
  integer, target :: i(10)
  integer, pointer :: j(:)
  i = [(n,n=1,10)]
  j => i
  print '(10(I3))', j
end program foo

When compiling with gfortran, there are no warnings/errors (there shouldn't be for a correct program), and the output is 1 2 3 4 5 6 7 8 9 10. A slightly more complicated version of the above is

program foo
  integer, pointer :: j(:)
  allocate(j(10))
  j = [(n,n=1,10)]
  print '(10(I3))', j
end program foo

The allocate statement will allocate what can be called an anonymous target with 10 elements (for the lack of a better name). j points at that anonymous target. The difference here compared to the first program is that j = [(n,n=1,10] is an intrinsic assignment while in the former j => i is pointer assignment.

PS: If one does not need a pointer, it is (IMO) good practice to use allocatable.