c - 调用 C_F_POINTER 时 Fortran 指针和 Fortran 可分配的区别

问题描述

问题是，“C_F_POINTER”以“可分配数组”作为参数成功编译（ifort 版本 19.0.5.281），并且在使用“指针”作为参数的情况下，它的工作方式完全相同。

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

在上面的测试代码中，我在代码的左侧将一些点标记为（*1）、（*2）等。

(*1) & (*2) : 代码的可分配数组版本和指针版本 (*3) : 我调用 'C_F_POINTER' (*4) : 打印语句以查看输出

在 (*1) 和 (*2) 情况下，结果“完全相同”，并且所有值都正确转换为所需的 CSR 格式。

 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

2 年前我在 StackOverflow 中发现了一个类似的问题（fortran 指针或 c_f_pointer 调用的可分配数组之间的区别）。

这个问题现在在我脑海中提出完全相同的问题。

如果我用我的话重新排列问题，

指针和可分配数组的区别？

据我所知，在 C 中，数组存储在连续的内存中，并且可以由指向其第一个元素的指针表示。在 Fortran90 中，如果我将数组作为“假定大小的数组”传递给子例程，代码的行为就像它从不关心它的分配方式、大小如何，并将数组视为存储在连续站点中的一维数组。
在下面的代码中，子程序 'assign_A' 只是将 'tot_array(1,2)' 作为其起点，并在连续的站点上执行它的工作，并且似乎甚至超出了 'tot_array' 的范围！（tot_array 是 2x2 矩阵，assign_A 的 do 循环从 tot_array(1,2) 开始运行 5 次）我“感觉”指针和可分配数组在这个意义上是类似的东西。但显然，作为c_f_pointer call 的 fortran 指针或可分配数组之间差异的答案，它们是不同的东西。为什么当数组作为“假定大小”传递给子例程时，它们就像指针一样？

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

在 Fortran90 中调用“C_F_POINTER”时使用“可分配数组”和“指针”有什么区别吗？
- 我使用了 ifort 版本 19.0.5.281，据我检查，这个编译器似乎给了我完全相同的结果。如果可以的话，我更喜欢使用可分配数组而不是指针。将“可分配数组”和“指针”与“C_F_POINTER”一起使用有什么区别，这样做有什么我应该注意的吗？
- 用于 c_f_pointer 调用的 fortran 指针或可分配数组之间差异的答案说我应该使用指针，而不是使用带有 C_F_POINTER 的可分配数组，但似乎这是一些持续存在的问题，当时还没有完全确定。为什么为 fortran 指针设计的“C_F_POINTER”可以与可分配数组一起正常工作并且结果相同，是否有任何结论？

感谢您阅读这个问题。

标签： cpointersfortranfortran-iso-c-binding

解决方案

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.

c - 调用 C_F_POINTER 时 Fortran 指针和 Fortran 可分配的区别

问题描述

解决方案

推荐阅读