argasm
argasm

Reputation: 379

Pointer array to an incontiguous array section

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?

Upvotes: 3

Views: 183

Answers (2)

argasm
argasm

Reputation: 379

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.

Upvotes: 0

IanH
IanH

Reputation: 21431

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.

Upvotes: 5

Related Questions