pointers, arrays of different dimensions, and derived types

Stephen_L_ · ‎12-27-2016

Hello! Right now I'm hoping to point data that is inside a derived type to an array (target) that is of different dimensions and is passed to the rest of the world (e.g. Python) via iso_c_binding. The application is a ray trace program, where each derived type is passed to and computed in a separate thread via openMP, and all the results are aggregated into the c type arrays. I'd like to point 1D arrays inside the derived types to subsets of higher-D (more than 3D) arrays.

The following is drawn from my current attempt to accomplish this. Obviously it doesn't work, so I'm grateful in advance either for tips to fix, or a possible alternate approach! (I feel like I'm headed the wrong way, though it's possibly because the shape that this code has taken over the last few years (my fault) has penned me in.)

module dataInterface

type rayType 
    real, dimension(:), pointer :: x, y, z
    real :: initElev, initAz, finElev, finAz, freq
end type

type inputParams 
    integer :: nPts, nFreq, nAz, nElev
end type 

end module

subroutine rayMain(inputParams, outputX, outputY, outputZ) bind(c, name="rayMain") 

use dataInterface
use iso_c_binding

implicit none

real(c_double), dimension(inputParams%nPts, inputParams%nFreq, inputParams%nElev, inputParams%nAz), target :: outputX, outputY, outputZ
type(rayType), dimension( inputParams%nFreq, inputParams%nElev, inputParams%nAz) :: rays
real(dp), dimension(inputParams%nFreq) :: f
real(dp), dimension(inputParams%nAz) :: az
real(dp), dimension(inputParams%nElev) :: elev
integer :: ii, jj, kk

do ii = 1,inputParams%nFreq
    do jj = 1, inputParams%nElev
        do kk = 1,inputParams%nAz 
            rays(ii,jj,kk)%x => outputX(:,ii,jj,kk)
            rays(ii,jj,kk)%y => outputY(:,ii,jj,kk)
            rays(ii,jj,kk)%z => outputZ(:,ii,jj,kk)
        enddo
    enddo
enddo

! this is not complete, obviously, but the ray trace is looped over the lists of frequencies and initial angles ...
call raytrace(rays, elev, az, f)

end subroutine

Note that there are other moving parts, such as porting onto GPGPU, so moving data onto and off of devices is a factor, as is overall memory usage efficiency (i.e., not carving off and writing/reading multiple copies of the same data in memory). I'm trying to narrow the scope of this question!

Stephen_L_ · ‎12-27-2016

UPDATE:

I got the above (forgive the minor errors, e.g. in kind type omissions) to build (ifort) with the additional changes:

type rayout

real(dp), dimension(:,:,:,:), pointer :: x, y, z

end type

rays%x(ii,jj,kk) => xout( : ,ii:ii,jj:jj,kk:kk)

etc.

But when I test it's throwing an access violation (0xFFFFFFFFFFFF) error.

In addition to possible fixes to this approach, I'd appreciate suggestions of another approach -- inasmuch as one can formulate such views based on the limited info presented here. Thanks!

Steven_L_Intel1 · ‎12-28-2016

A derived type with POINTER components is not interoperable with other languages. ISO_C_BINDING is just a set of declarations, not a method. If you want to have the compiler help you with interoperability, use BIND(C) on your structures and interfaces. Perhaps what you need here is to have the x, y and z components be of type C_PTR and assign them the C_LOC() of your data. Then the other languages will just see a data pointer.

FortranFan · ‎12-29-2016

Stephen L. wrote:

.. But when I test it's throwing an access violation (0xFFFFFFFFFFFF) error.

In addition to possible fixes to this approach, I'd appreciate suggestions of another approach -- inasmuch as one can formulate such views based on the limited info presented here. Thanks!

@Stephen L.,

You don't provide any reproducer of your issue(s), so it's anyone's guess what the exact problem(s) might be. However, looking at your original post and the subsequent comment, I wonder if you are facing several hurdles:

1) with gaining a full understanding of the facilities available in current Fortran standard, as supported by Intel Fortran compiler; toward this aspect, this forum and Dr Fortran articles and especially the references in this link might be of help,

2) with differences in Fortran and with the other languages that will call your procedures; again, some research into Python and C/C++ might be due, especially with row-major versus column-major array storage mechanisms, etc., and

3) with program design and setting up of Fortran library code with appropriate APIs and so forth; here also, more investigation will help but what might be even more important is to setup simple prototypes that enact your actual code design needs in small chunks; you can then work your way up to your full-blown Fortran library design.

It would appear the main need expressed in your original post is that some program with a companion C processor will call your Fortran code, the caller will have 'data' allocated in multi-dimensional arrays, and you want the Fortran to make use of such 'data' without having to allocate additional memory in the Fortran library, and toward this you want to make use of POINTER facility in Fortran. All well and good, except if you are dealing with multidimensional arrays then you need to pay attention to storage sequences on both sides of the interface and if you are making use of POINTERs, you need to be very careful with dangling pointers and unreachable memory locations and so forth. You don't show how you take care of such details and any gaps in these aspects can cause major problems including access violations.

My suggestion will be for you to avoid POINTERs if you can; but if you must, then make use of OO concepts to further extend your data types (e.g., rayType) and utilize the concept of information hiding to only provide controlled routes (setter methods) to setup data pointers, and further implement finalizer(s) to clean up pointer references, etc. Then setup small working examples and try them out until they function as intended, ensure there are no memory leaks (Intel Inspector can assist) and work your way up. Shown below is a go-by, any of the references in the above link can help explain the language details if they are unclear:

module data_m

   use, intrinsic :: iso_c_binding, only : c_int

   implicit none

   private

   type, public :: data_t
      private
      integer(c_int), pointer :: m_dat1(:) => null()
      integer(c_int), pointer :: m_dat2(:) => null()
   contains
      private
      final :: clean_data_t
      procedure, pass(this) :: set_data_t
      procedure, pass(this) :: calc_data_t
      generic, public :: set => set_data_t
      generic, public :: calc => calc_data_t
   end type data_t

contains

   subroutine clean_data_t( this )

      !.. Argument list
      type(data_t), intent(inout) :: this

      !.. Perform clean-up operations
      if ( associated(this%m_dat1) ) then
         this%m_dat1 => null()
      end if
      if ( associated(this%m_dat2) ) then
         this%m_dat2 => null()
      end if

      return

   end subroutine clean_data_t

   subroutine set_data_t( this, dat )

      !.. Argument list
      class(data_t), intent(inout)          :: this
      integer(c_int), intent(inout), target :: dat(:,:)

      !.. Set up data pointers
      this%m_dat1 => dat(:,1)
      this%m_dat2 => dat(:,2)

      return

   end subroutine set_data_t

   subroutine calc_data_t( this )

      !.. Argument list
      class(data_t), intent(inout) :: this

      !.. Perform some calcs
      print *, "In data_m%calc_data_t:"
      print *, "upon entry: this%m_dat1 = ", this%m_dat1
      print *, "            this%m_dat2 = ", this%m_dat2

      this%m_dat1 = -this%m_dat1
      this%m_dat2 = -this%m_dat2

      print *, "Following calculations: this%m_dat1 = ", this%m_dat1
      print *, "                        this%m_dat2 = ", this%m_dat2

      return

   end subroutine calc_data_t

end module data_m

module api_m

   use, intrinsic :: iso_c_binding, only : c_int, c_size_t, c_ptr, c_f_pointer
   use data_m, only : data_t

   implicit none

   private

   type(data_t), allocatable, save :: gdat  ! reusable variable to work space variable

   public :: main_api
   public :: cleanup_api

contains

   subroutine main_api( dim1, dim2, dat_ptr, irc ) bind(C, name="main_api")

      !.. Argument list
      integer(c_size_t), intent(in), value :: dim1
      integer(c_size_t), intent(in), value :: dim2
      type(c_ptr), intent(in), value       :: dat_ptr
      integer(c_int), intent(inout)        :: irc

      irc = 0
      if ( dim2 < 2 ) then
         irc = 1
         ! Elided is the check for appropriate dimensions
         return
      end if

      !.. Allocate work space variable if needed
      if ( .not. allocated(gdat) ) then
         allocate( gdat, stat=irc )
         if ( irc /= 0 ) then
            !.. error handling elided
            return
         end if
      end if

      !.. Setup data pointers
      blk: block
         integer(c_int), pointer :: dat(:,:) => null()
         call c_f_pointer( cptr=dat_ptr, fptr=dat, shape=[ dim1, dim2 ] )
         call gdat%set( dat )
         dat => null()
      end block blk

      !.. Perform calcs
      call gdat%calc()

      return

   end subroutine main_api

   subroutine cleanup_api( irc ) bind(C, name="cleanup_api")

      !.. Argument list
      integer(c_int), intent(inout) :: irc

      irc = 0
      !.. Deallocate data if needed
      if ( allocated(gdat) ) then
         deallocate( gdat, stat=irc )
         if ( irc /= 0 ) then
            !.. error handling elided
            return
         end if
      end if

      return

   end subroutine cleanup_api

end module api_m

#include <stdio.h>

extern void main_api(size_t, size_t, int *, int *);
extern void cleanup_api(int *);

int main()
{
   enum PROB_SIZE { N = 3, M = 2 };

   int dat;
   int i;
   int j;
   int irc;

   for (i = 0; i < M; i++) {
      for (j = 0; j < N; j++) {
         dat = i*N + j + 1;
      }
   }
   printf("In caller, initially:\n");
   for (i = 0; i < M; i++) {
      for (j = 0; j < N; j++) {
         printf("dat[%d][%d] = %d\n", i, j, dat);
      }
   }

   main_api(N, M, &dat[0][0], &irc);
   if (irc != 0) {
      printf("main_api returned an error: irc = %d", irc);
      return(irc);
   }

   cleanup_api(&irc);
   if (irc != 0) {
      printf("cleanup_api returned an error: irc = %d", irc);
      return(irc);
   }

   printf("In caller, following Fortran invocation:\n");
   for (i = 0; i < M; i++) {
      for (j = 0; j < N; j++) {
         printf("dat[%d][%d] = %d\n", i, j, dat);
      }
   }

   return 0;
}

Upon execution with Intel Fortran,

In caller, initially:
dat[0][0] = 1
dat[0][1] = 2
dat[0][2] = 3
dat[1][0] = 4
dat[1][1] = 5
dat[1][2] = 6
 In data_m%calc_data_t:
 upon entry: this%m_dat1 =  1 2 3
             this%m_dat2 =  4 5 6
 Following calculations: this%m_dat1 =  -1 -2 -3
                         this%m_dat2 =  -4 -5 -6
In caller, following Fortran invocation:
dat[0][0] = -1
dat[0][1] = -2
dat[0][2] = -3
dat[1][0] = -4
dat[1][1] = -5
dat[1][2] = -6
Press any key to continue . . .