fortran 2003 and pointers

Question

I am writing a Fortran program. The program implements some numerical methods. Program speed is very important. I decided to get rid of dynamic arrays (whether it speeds up the program?), and faced with the following problem. I have 3d arrays (NXxNYxNZ = MAX elements) where I know MAX, but I dont know NX/NY/NZ ratio. It can be like this : 1x1xNZ or like this 2xNYx1 , etc. The solution I see - use pointers. Simplified 2D case:

program ptrtest
   parameter ( MAX = 50 )    ! I know this number on compile time.
   integer :: NX,NY          ! I will find this numbers on run time.
   real, target :: a(MAX)    ! Static Array
   real, pointer :: b(:,:)
   a = 5
   read(*,*) NX, NY          ! I get NX, NY (NX*NY == MAX)
   b (1:NX, 1:NY) => a       ! I can use b(:,:) <- this is my goal.
end program ptrtest

This example works, but i am afraid that such update slow down the real program where I am using even 5d arrays. Is it possible?

Answer 1

You say that program speed is important, so old-style Fortran will give you the best:

• allocate a 1D array of size MAX
• pass is to a subroutine as well as NX,NY,and NZ. This will avoid the pointer indirection that will give you a loss of performance (unlike C, the Fortran standard assumes that subroutine array arguments don't overlap, see this post).

For example:

---

program noptrtest
   implicit none
   integer, parameter :: MAX = 50     ! I know this number on compile time.
   integer :: NX,NY          ! I will find this numbers on run time.
   real    :: a(MAX)         ! Static Array

   read(*,*) NX, NY          ! I get NX, NY (NX*NY == MAX)

   if (NX*NY /= MAX) then
     stop 'NX*NY /= MAX'
   endif
   call run(a,NX,NY)
end program noptrtest


subroutine run(b,NX,NY)
  integer :: NX,NY
  real    :: b(NX,NY)
  integer :: i,j

  do j=1,NY
    do i=1,NX
      ! use b(i,j) here
    enddo
  enddo
end

---

If performance really matters here are some other useful tricks:

• Tell your compiler to align the arrays on a 32-byte boundary (with ifort, use the !DIR$ ATTRIBUTES ALIGN : 32 directive
• Dimension physically your 2D array such that the size of leading dimension is a multiple of 32-bytes. For real*4, you need a multiple of 8 elements.
• Tell the compiler that each column is properly aligned using the !DIR$ VECTOR ALIGNED directive

---

For example:

program noptrtest
   implicit none
   integer, parameter :: real_kind = 4 
   integer, parameter :: MAX = 50               ! I know this number on compile time.
   integer, parameter :: MAXDIM = MAX*(32/real_kind)    ! Max possible dimension required
   integer            :: NX,NY, NX_MOD         ! I will find this numbers on run time.
   real(real_kind)    :: a(MAXDIM)             ! Static Array

   !DIR$ ATTRIBUTES ALIGN : 32 :: a

   read(*,*) NX, NY          ! I get NX, NY (NX*NY == MAX)

   if (NX*NY /= MAX) then
     stop 'NX*NY /= MAX'
   endif

   if (mod(NX,real_kind) == 0) then
     NX_MOD = NX
   else
     NX_MOD = ((NX/real_kind)+1)*real_kind
   endif

   call run(a,NX_MOD,NX,NY)
end program noptrtest


subroutine run(b,NX_MOD,NX,NY)
  integer :: NX_MOD,NX,NY
  real    :: b(NX_MOD,NY)
  integer :: i,j

  do j=1,NY
    !DIR$ VECTOR ALIGNED
    do i=1,NX
    ! use b(i,j) here
    enddo
  enddo
end

---

Edit

This old-style Fortran trick avoids pointer aliasing.

References on pointer aliasing:

• http://en.wikipedia.org/wiki/Aliasing_(computing) , section "Conflicts with optimization"
• gcc C/C++ assume no pointer aliasing
• Performance impact of -fno-strict-aliasing
• http://scc.qibebt.cas.cn/docs/compiler/intel/11.1/Intel%20C++%20Compiler%20for%20Linux/main_cls/optaps/common/optaps_perf_run.htm