Difference in Fortran pointer and Fortran allocatable in calling C_F_POINTER

Question

The thing is, 'C_F_POINTER' compiles successfully(ifort version 19.0.5.281) with 'allocatable arrays' as its argument, and it works in the exactly same way with the case in which 'pointer' is used as its argument.

program test1
    use mkl_spblas
    use omp_lib
    use iso_c_binding

    implicit none
    integer, parameter    :: DIM_ = 4, DIM_2 = 6
    integer               :: stat, i
    integer               :: irn(DIM_2), jcn(DIM_2)
    real*8                :: val(DIM_2)
    integer(c_int)        :: indexing
    integer               :: DIM_r, DIM_c
    type(c_ptr)           :: rows_start_c, rows_end_c, col_indx_c, values_c
(*1)!integer,allocatable   :: rows_start_f(:), rows_end_f(:), col_indx_f(:)
    !real*8 ,allocatable   :: values_f(:)
(*2)integer  ,pointer     :: rows_start_f(:), rows_end_f(:), col_indx_f(:)
    real*8   ,pointer     :: values_f(:)
    type(SPARSE_MATRIX_T) :: mat1, mat2

    irn = (/ 2, 2, 3, 4, 0, 0 /)
    jcn = (/ 1, 2, 3, 2, 0, 0 /)
    val = (/ 5, 8, 3, 6, 0, 0 /)

    call omp_set_num_threads(1)

    stat = mkl_sparse_d_create_coo (A=mat1, indexing=SPARSE_INDEX_BASE_ONE, &
                                    rows=DIM_, cols=DIM_, nnz=DIM_,&
                                    row_indx=irn, col_indx=jcn, values=val  )
    if (stat /= 0) stop 'Error in mkl_sparse_d_create_coo'

    stat = mkl_sparse_convert_csr (source=mat1,&
                                   operation=SPARSE_OPERATION_NON_TRANSPOSE, &
                                   dest = mat2 )
    if (stat /= 0) stop 'Error in mkl_sparse_convert_csr'

    stat = mkl_sparse_d_export_csr(mat2, indexing, DIM_r, DIM_c,  &
                                   rows_start_c, rows_end_c, col_indx_c, values_c)

(*3)call c_f_pointer(rows_start_c, rows_start_f, [DIM_r])
    call c_f_pointer(rows_end_c  , rows_end_f  , [DIM_c])
    call c_f_pointer(col_indx_c  , col_indx_f  , [rows_end_f(DIM_r)-1])
    call c_f_pointer(values_c    , values_f    , [rows_end_f(DIM_r)-1])

    stat = mkl_sparse_destroy (A=mat1)
    if (stat /= 0) stop 'Error in mkl_sparse_destroy (mat1)'

    stat = mkl_sparse_destroy (A=mat2)
    if (stat /= 0) stop 'Error in mkl_sparse_destroy (mat2)'

    call mkl_free_buffers

(*4)print *, 'rows_start'
    print *, rows_start_f
    print *, 'rows_end'
    print *, rows_end_f
    print *, 'col_indx'
    print *, col_indx_f
    print *, 'values'
    print *, values_f
    print *, 'indexing'
    print *, indexing
    print *, 'size(values_f,1)'
    print *, size(values_f,1)

end program test1

In the test code above, I marked some points as (*1), (*2), and so on in the leftside of the code.

(*1) & (*2) : allocatable array version and pointer version of the code (*3) : where I call 'C_F_POINTER' (*4) : print statements to see the output

The results are 'exactly' the same in both (*1), and (*2) case, and all values are properly converted into desired CSR format.

 rows_start
           1           1           3           4
 rows_end
           1           3           4           5
 col_indx
           1           2           3           2
 values
   5.00000000000000        8.00000000000000        3.00000000000000     
   6.00000000000000     
 indexing
           1
 size(values_f,1)
           4

I found a similar question in StackOverflow 2 years ago (difference between fortran pointers or allocatable arrays for c_f_pointer call).

This question is asking the exactly the same questions in my mind right now.

If I rearange questions in my words,

Difference between pointer and allocatable array?

In C, as far as I know, the arrays are stored in contiguous memory and can be represented by the pointer which points its 1st element. And in Fortran90, if I pass a array into a subroutine as 'assumed-size array', the code behaves like it never cares about how it's allocated, how it's size is like, and treates the array as 1D being stored in contiguous site.
In below code, the subroutine 'assign_A' just gets the 'tot_array(1,2)' as its starting point, and do its work on contiguous site and seems to do it even out of bound of 'tot_array'!! (tot_array is 2x2 matrix, and assign_A's do loop runs 5 times starting at tot_array(1,2)) I was 'feeling' the pointer and allocatable arrays are similar stuff in this sense. But apparently, as the answers in difference between fortran pointers or allocatable arrays for c_f_pointer call, they are different things. Why arrays acts like pointer when they are passed to subroutine as 'assumed-size' one?

program assumed_size_array_test
  implicit none
  external assign_A
  real*8 :: tot_array(2,2)
  integer:: i
  
  ! Initially 'tot_array' set to be 1.d0
  tot_array = 1.d0
  
  write(*,*) 'Before'
  write(*,'(5f5.2)') tot_array
  
  call assign_A(tot_array(1,2))
  
  write(*,*) 'After'
  write(*,'(5f5.2)') tot_array

end program

subroutine assign_A(A)
  implicit none
  real*8, intent(inout) :: A(*)
  integer :: i
  
  do i = 1,5
    A(i) = 2.d0
  enddo

end subroutine

 Before
 1.00 1.00 1.00 1.00
 After
 1.00 1.00 2.00 2.00

Is there any difference in using 'allocatable array' and 'pointer' in calling 'C_F_POINTER' in Fortran90?
- I used ifort version 19.0.5.281, and this compiler seems to give me exactly the same results as far as I checked. If it's okay, I prefer to use allocatble arrays instead of pointers. Is there any difference in using 'allocatable array' and 'pointer' with 'C_F_POINTER', and is there anything that I should be aware of in doing so?
- The answers in difference between fortran pointers or allocatable arrays for c_f_pointer call says that I SHOULD use pointers, not using allocatable arrays with C_F_POINTER, but it seems it's some ongoing issue that was not concluded exactly at that time. Is there any conclusion in why 'C_F_POINTER', which is designed for fortran pointer, works fine with allocatable arrays and is result is the same?

Thank you for reading this question.

The comments in the linked question appear to apply here also: if you give c_f_pointer an allocatable rather than a pointer, what you have isn't a Fortran program. (And the compiler can do whatever it likes with it.) — francescalus
– francescalus, Commented May 13, 2020 at 14:48

Vladimir F Героям слава · Accepted Answer · 2020-05-13 15:28:55Z

Obviously, both Fortran POINTER variables and ALLOCATABLE variables have a lot of common in their internal impementation. Most of that is under the hood and should not be accessed directly. Both allocate some memory and probably use the same operating system's or C runtime library's allocator. For example, malloc().

In both there is some memory allocated or pointed to and described by a simple address (for scalars) or by an array descriptor (for an array).

Pointers and allocatable variables mainly differ in what you can do with them and what the compiler will do with them for you. You can think of allocatables as a sort of "smart pointers" quite similar to std::unique_ptr in C++. Recall what happens in C++ you have new and delete which in turn call malloc and free but you are not allowed to mix them. And you are certainly not allowed to manually modify the address stored in a C++ smart pointer either.

When you send an allocatable variable to a procedure that expects a pointer, anything can happen, it is an undefined behaviour. But, if the internal hidden structure has a similar layout, it may happen that you actually set the allocatable internals to point to some memory that was not allocated through allocatable. You may then think that everything is OK and you have a new feature. However, when the time for deallocation comes, and allocatables are often deallocated automatically, it can easilly fail in very unpredictable ways. It can crash in very strange places of the code, the results can be wrong and so on. Anything can happen.

For example, this extremely ugly program works for me too (in gfortran):

subroutine point(ptr, x)
  pointer :: ptr
  target :: x
  ptr => x
end subroutine

  interface
    subroutine point(ptr, x)
      allocatable :: ptr
      target :: x
    end subroutine    
  end interface

  allocatable z

  y = 1.0

  call point(z, y)

  print *, z
end

But you should never do stuff like this. It is really something very, very wrong. If you make z a local variable, so that it is deallocated, or if you try to deallocate it, it will crash. That is because the only information the compiler has is the address. Internally, the allocatable really looks the same as a pointer. It is just an address (for a scalar). The only difference is what you are allowed to do with it and what the compiler will do for you automatically.

This won't even crash, because the internal implementation similarities I mentioned. but it is no less wrong.

subroutine point(ptr, x)
  pointer :: ptr
  target :: x
  ptr => x
end subroutine     

  interface
    subroutine point(ptr, x)
      allocatable :: ptr
      target :: x
    end subroutine    
  end interface

  allocatable z
  pointer y

  allocate(y)
  y = 1.0

  call point(z, y)

  print *, z

  deallocate(z)
end

It just survives because both allocatable and pointer use the same internal allocator (malloc) in gfortran and they are both implemented as a simple address.

Thank you for such a detailed answer! It helps a lot to understand how compiler, pointer, and allocatables do their work :-) Have a nice day

Collectives™ on Stack Overflow

Difference in Fortran pointer and Fortran allocatable in calling C_F_POINTER

1 Answer 1

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related