fftw3: Allocating aligned memory in Fortran

 
 7.5 Allocating aligned memory in Fortran
 ========================================
 
 In order to obtain maximum performance in FFTW, you should store your
 data in arrays that have been specially aligned in memory (SeeSIMD
 alignment and fftw_malloc).  Enforcing alignment also permits you to
 Functions::) to apply a given plan to more than one pair of in/out
 arrays.  Unfortunately, standard Fortran arrays do _not_ provide any
 alignment guarantees.  The _only_ way to allocate aligned memory in
 standard Fortran is to allocate it with an external C function, like the
 'fftw_alloc_real' and 'fftw_alloc_complex' functions.  Fortunately,
 Fortran 2003 provides a simple way to associate such allocated memory
 with a standard Fortran array pointer that you can then use normally.
 
    We therefore recommend allocating all your input/output arrays using
 the following technique:
 
   1. Declare a 'pointer', 'arr', to your array of the desired type and
      dimensions.  For example, 'real(C_DOUBLE), pointer :: a(:,:)' for a
      2d real array, or 'complex(C_DOUBLE_COMPLEX), pointer :: a(:,:,:)'
      for a 3d complex array.
 
   2. The number of elements to allocate must be an 'integer(C_SIZE_T)'.
      You can either declare a variable of this type, e.g.
      'integer(C_SIZE_T) :: sz', to store the number of elements to
      allocate, or you can use the 'int(..., C_SIZE_T)' intrinsic
      function.  e.g.  set 'sz = L * M * N' or use 'int(L * M * N,
      C_SIZE_T)' for an L x M x N array.
 
   3. Declare a 'type(C_PTR) :: p' to hold the return value from FFTW's
      allocation routine.  Set 'p = fftw_alloc_real(sz)' for a real
      array, or 'p = fftw_alloc_complex(sz)' for a complex array.
 
   4. Associate your pointer 'arr' with the allocated memory 'p' using
      the standard 'c_f_pointer' subroutine: 'call c_f_pointer(p, arr,
      [...dimensions...])', where '[...dimensions...])' are an array of
      the dimensions of the array (in the usual Fortran order).  e.g.
      'call c_f_pointer(p, arr, [L,M,N])' for an L x M x N array.
      (Alternatively, you can omit the dimensions argument if you
      specified the shape explicitly when declaring 'arr'.)  You can now
      use 'arr' as a usual multidimensional array.
 
   5. When you are done using the array, deallocate the memory by 'call
      fftw_free(p)' on 'p'.
 
    For example, here is how we would allocate an L x M 2d real array:
 
        real(C_DOUBLE), pointer :: arr(:,:)
        type(C_PTR) :: p
        p = fftw_alloc_real(int(L * M, C_SIZE_T))
        call c_f_pointer(p, arr, [L,M])
        _...use arr and arr(i,j) as usual..._
        call fftw_free(p)
 
    and here is an L x M x N 3d complex array:
 
        complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:)
        type(C_PTR) :: p
        p = fftw_alloc_complex(int(L * M * N, C_SIZE_T))
        call c_f_pointer(p, arr, [L,M,N])
        _...use arr and arr(i,j,k) as usual..._
        call fftw_free(p)
 
    See SeeReversing array dimensions for an example allocating a
 single array and associating both real and complex array pointers with
 it, for in-place real-to-complex transforms.