Pointer array to an incontiguous array section

Question

Can I have a pointer array that points to an incontiguous section of a target array? In the following program:

program test
  implicit none

  integer :: i
  integer, dimension(4), target :: orig
  integer, dimension(:), pointer :: ref

  orig = (/ (i , i = 1,4) /)

  ref(1:2) => orig(1:2)   
  ref(3:3) => orig(4:4)

  print *, 'orig = ', orig
  print *, 'ref = ', ref

end program test

I want the pointer ref to point to the incontiguous section (/1,2,4/) of the target array orig. The code compiles without a warning, however the output is not desired:

code output:
 orig =            1           2           3           4
 ref =            4

It appears as the compiler forces the pointer to point to a contiguous subsection even though I first point to indices 1:2 but compiler only forces pointing to 4:4. As a result my pointer at the end only points to one entry.

I can pack the entries (/1,2,4/) into another array and then use it contiguously but this would not be efficient Fortran coding for large data set. Do you have any idea how to make this working?

IanH · Accepted Answer · 2014-08-07 22:59:00Z

5

Your pointer assignments have what is known as bounds remapping. The "entire" pointer object is pointer assigned, but it is given the bounds indicated in the parentheses on the left hand side. This lets you have a lower bound other than one and have the rank of the pointer differ from the target. This feature is similar to the capabilities that you have with array argument association when the dummy array is explicit or assumed size.

You cannot have a pointer point to arbitrary elements of a target array - the elements in the target have to be "regular" - have a constant stride between them.

As an alternative to repacking you can obviously store those indices in an array and use the combination of the whole target and that array of indices to get the behaviour that you want, perhaps by wrapping a pointer to the whole target and the array of indices in a derived type with defined assignments.

answered Aug 7, 2014 at 22:59

IanH

21.6k2 gold badges41 silver badges61 bronze badges

Sign up to request clarification or add additional context in comments.

6 Comments

argasm Over a year ago

In C this is very easy to do with a double pointer. Regarding your suggestion, I'm not sure if I define my own wrapping and overloaded operators I will get the same performance like the native array operations in Fortran give me. Pointing to nonuniform sections seems to be very easy to be implemented in a Fortran compiler. I wonder why compiler developers hasn't done such thing yet. This is the first time I feel "C" is better for number crunching!

IanH Over a year ago

You are confusing simplicity of source code with performance.

Vladimir F Героям слава Over a year ago

The compiler developers implement the Fortran standard, not their arbitrary thougts what might be useful. In your situation you may try vector indexing. The packed array could be actually pretty efficient for vectorization.

argasm Over a year ago

@VladimirF Thanks! What do you mean by vector indexing? Is it like orig((/1, 2, 4/)) in this example?

argasm Over a year ago

@VladimirF so I guess if I use vector indexing like you said, for example orig((/1, 2, 4/)) would pack the entries in a contiguous temporary array and this is not what I want. I can simply do this kind of vector indexing manually with one for loop and the reason I never use this kind of indexing is I am worried about temporary packed copies. Correct me if I'm wrong. Thanks

|

argasm · Accepted Answer · 2014-08-08 19:57:30Z

Following "@IanH" suggestion, I implemented a wrapper which only has pointers to indices and data and everything is done IN PLACE so this should be expected to be more appropriate for large data. For small data set, I think "@Vladimir F" idea of vector indexing is better and more readable. The following code passed valgrind full memleak test and I tried for large data set (6GB double precision) and it doesn't copy data.

module any_stride
  implicit none
  private
  ! ------------------------------------------------------------
  ! data = pointer to the contiguous 1D double precision array
  ! indx = pointer to the non-uniform indices
  ! example: for a stride where: 
  ! stride%data => (/ 1.0d0, 2.0d0, 3.0d0, 4.0d0 /)
  ! stride%indx => (/ 1, 4 /)
  ! is a wrapper around 1d data array with four elements
  ! that only elements 1 and 4 are visible to the outside world  
  !
  type stride
     private
     real*8, dimension(:), pointer, public :: data => null()
     integer, dimension(:), pointer, public :: indx => null()
   contains
     procedure :: assgn_from_scalar
     procedure :: assgn_from_array
     ! ... 
     ! ...
     ! ... add more overloaded procs if you want 
     generic :: assignment(=) => assgn_from_scalar &
          , assgn_from_array
     generic :: operator(+) => add_to_scalar &
          , add_to_array
     ! ... 
     ! ...
     ! ... add more overloaded procs if you want 
  end type stride

  public :: stride

contains

   subroutine assgn_from_scalar(this, rhs)
    implicit none
    class(stride), intent(inout) :: this
    real*8, intent(in) :: rhs

    ! local vars
    integer :: i

    do i = 1, size(this%indx)
       this%data(this%indx(i)) = rhs
    end do

    ! done here
  end subroutine assgn_from_scalar

  subroutine assgn_from_array(this, rhs)
    implicit none
    class(stride), intent(inout) :: this
    real*8, dimension(:), intent(in) :: rhs

    ! local vars
    integer :: i

    do i = 1, size(this%indx)
       this%data(this%indx(i)) = rhs(i)
    end do

    ! done here
  end subroutine assgn_from_array


  function add_to_scalar(lhs, rhs)
    implicit none
    class(stride), intent(in), target :: lhs
    real*8, intent(in) :: rhs
    class(stride), pointer :: add_to_scalar

    ! local vars
    integer :: i

    do i = 1, size(lhs%indx)
       lhs%data(lhs%indx(i)) = lhs%data(lhs%indx(i)) + rhs
    end do

    add_to_scalar => lhs

    ! done here
  end function add_to_scalar

  function add_to_array(lhs, rhs)
    implicit none
    class(stride), intent(in), target :: lhs
    real*8, dimension(:), intent(in) :: rhs
    class(stride), pointer :: add_to_array

    ! local vars
    integer :: i

    do i = 1, size(lhs%indx)
       lhs%data(lhs%indx(i)) = lhs%data(lhs%indx(i)) + rhs(i)
    end do

    add_to_array => lhs

    ! done here
  end function add_to_array

end module any_stride

! a little tester program
!
! COMMENT and UNCOMMENT after this line 
! when using the module in this file
!
program test
  use any_stride
  implicit none

  ! local vars
  integer :: i
  real*8, dimension(4), target :: data
  real*8, dimension(2) :: rhs = (/ 3.0d0, 6.0d0 /)
  integer, dimension(2), target :: indx
  type(stride) :: astride

  data = (/ 1.0d0, 2.0d0, 3.0d0, 4.0d0 /)
  indx = (/ 1, 4 /)
  astride = stride(data, indx)


  astride = astride + rhs
  print *, 'astride data= ', astride%data

  ! done here
end program test

the values 1 and 4 are only added to the first and the last entries of array data.

Collectives™ on Stack Overflow

Pointer array to an incontiguous array section

2 Answers 2

6 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

6 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related