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

Question

事实是，“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 中发现了一个类似的问题 ( difference between fortran pointers or allocatable arrays for c_f_pointer call )。

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

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

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

• 据我所知，在 C 中，数组存储在连续的内存中，可以用指向其第一个元素的指针来表示。在 Fortran90 中，如果我将一个数组作为“假定大小的数组”传递给子例程，代码的行为就像它从不关心它是如何分配的，它的大小是怎样的，并将数组视为存储在连续站点中的一维数组。
• 在下面的代码中，子例程“assign_A”只是将“tot_array(1,2)”作为其起点，并在连续的站点上进行工作，而且似乎甚至超出了“tot_array”的范围!! (tot_array 是 2x2 矩阵，assign_A 的 do 循环从 tot_array(1,2) 开始运行 5 次)我“感觉”指针和可分配数组在这个意义上是相似的东西。但显然，正如 difference between fortran pointers or allocatable arrays for c_f_pointer call 中的答案一样，它们是不同的东西。为什么数组在作为“假定大小”传递给子例程时表现得像指针？

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

1. 在 Fortran90 中调用“C_F_POINTER”时使用“可分配数组”和“指针”有什么区别吗？
   • 我使用的是 ifort 版本 19.0.5.281，据我检查，这个编译器似乎给出了完全相同的结果。如果可以的话，我更喜欢使用 allocatble 数组而不是指针。将“可分配数组”和“指针”与“C_F_POINTER”一起使用有什么区别吗？这样做有什么我应该注意的吗？
   • difference between fortran pointers or allocatable arrays for c_f_pointer call 中的答案说我应该使用指针，而不是使用带有 C_F_POINTER 的可分配数组，但它似乎是一些持续存在的问题，当时还没有完全确定。对于为什么为 Fortran 指针设计的“C_F_POINTER”在可分配数组上工作良好并且结果相同，是否有任何结论？

感谢您阅读这个问题。

Answer 1

显然，Fortran POINTER 变量和ALLOCATABLE 变量在内部实现上有很多共同点。其中大部分都在幕后，不应直接访问。两者都分配一些内存并且可能使用相同的操作系统或 C 运行时库的分配器。例如，malloc()。

两者都分配或指向一些内存，并通过简单地址(对于标量)或数组描述符(对于数组)进行描述。

指针和可分配变量的主要区别在于您可以使用它们做什么以及编译器将为您做什么。您可以将可分配项视为一种与 C++ 中的 std::unique_ptr 非常相似的“智能指针”。回想一下 C++ 中发生的事情，你有 new 和 delete，它们依次调用 malloc 和 free 但你不允许混合它们。而且您当然也不允许手动修改存储在 C++ 智能指针中的地址。

当您将可分配变量发送到需要指针的过程时，任何事情都可能发生，这是未定义的行为。但是，如果内部隐藏结构具有类似的布局，则可能会发生您实际上将 allocatable internals 设置为指向一些未通过 allocatable 分配的内存。然后您可能会认为一切正常并且您拥有了一项新功能。然而，当释放时间到来时，可分配对象通常会自动释放，它很容易以非常不可预测的方式失败。它可能会在代码非常奇怪的地方崩溃，结果可能是错误的等等。任何事情都有可能发生。

例如，这个极其丑陋的程序也适用于我(在 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

但是你不应该永远做这样的事情。这真的是非常非常错误的事情。如果您将 z 设为局部变量，以便它被释放，或者如果您尝试释放它，它将崩溃。那是因为编译器拥有的唯一信息是地址。在内部，可分配对象看起来真的和指针一样。它只是一个地址(对于标量)。唯一的区别是您可以用它做什么以及编译器会自动为您做什么。

这甚至不会崩溃，因为我提到的内部实现相似性。但这同样是错误的。

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

它幸存下来是因为可分配和指针在 gfortran 中使用相同的内部分配器 (malloc)，并且它们都作为简单地址实现。