Fortran中子程序调用传递切片数组时如何避免创建临时数组？

问题描述 #

我编写了如下Fortran子程序：

subroutine invertCell( a1vec,a2vec,a3vec, rlvec ) !(icluster)
    implicit none
    ! ======= Parameters
    real, dimension (3),  intent(in)  :: a1vec, a2vec, a3vec
    real, dimension (3,3),intent(out) :: rlvec
    ! ... Body ...
end subroutine

当我从3×3矩阵中按行切片作为参数调用该子程序时：

real, dimension (3,3) :: ervec
real, dimension (3,3) :: elvec
...
call invertCell( ervec(1,:),ervec(2,:),ervec(3,:), elvec )

会收到警告：Fortran runtime warning: An array temporary was created。

我可以通过显式创建临时数组规避警告，但写法很繁琐：

real, dimension (3)     :: ervec1,ervec2,ervec3
...
ervec1(:) = ervec(1,:)
ervec2(:) = ervec(2,:)
ervec3(:) = ervec(3,:)
call invertCell( ervec1,ervec2,ervec3, elvec )

我有C/C背景，C/C中数组可通过指针算术实现无复制切片，疑惑Fortran为何无法做到。

原因分析 #

Fortran采用列主序存储数组，即矩阵的列是连续内存块，而行是分散的非连续内存。你的子程序形参a1vec是固定大小的连续一维数组，当传入ervec(1,:)（矩阵的第一行，非连续内存）时，编译器为了满足形参对连续内存的要求，会自动创建临时数组，把切片数据复制进去，这就是警告的来源。

而C/C++是行主序存储，矩阵的行是连续内存，且函数参数可直接接受指针，因此能直接传递行切片的地址，无需复制。

无需显式临时数组的解决方案 #

1. 修改子程序为接受假定形状数组（Fortran 2003+） #

将形参改为假定形状数组（dimension(:)），这类数组的描述符包含内存布局信息，可直接处理非连续切片，无需创建临时数组：

subroutine invertCell( a1vec,a2vec,a3vec, rlvec ) !(icluster)
    implicit none
    ! ======= Parameters
    real, dimension(:),  intent(in)  :: a1vec, a2vec, a3vec
    real, dimension(3,3),intent(out) :: rlvec
    ! 可选：添加输入长度校验，避免非法参数
    if (size(a1vec)/=3 .or. size(a2vec)/=3 .or. size(a3vec)/=3) then
        error stop "invertCell: Input vectors must be length 3"
    end if
    ! ... Body ...
end subroutine

修改后直接用原调用方式即可，不会再产生临时数组警告。

2. 调整数组存储为列主序匹配形式 #

如果允许调整数据存储方式，将矩阵的行数据转为列存储，这样切片列时是连续内存，直接传入即可：

real, dimension(3,3) :: ervec, elvec
! 将原行数据转存为列（示例）
ervec(:,1) = [1.0, 2.0, 3.0]  ! 原第一行
ervec(:,2) = [4.0, 5.0, 6.0]  ! 原第二行
ervec(:,3) = [7.0, 8.0, 9.0]  ! 原第三行

! 调用时传入列切片（连续内存）
call invertCell( ervec(:,1), ervec(:,2), ervec(:,3), elvec )

3. 使用指针关联切片（Fortran 90+） #

用指针直接关联矩阵的行切片，指针会指向原数组内存，无需复制：

real, dimension(3,3) :: ervec, elvec
real, dimension(:), pointer :: er1, er2, er3

! 指针关联到矩阵的行
er1 => ervec(1,:)
er2 => ervec(2,:)
er3 => ervec(3,:)

call invertCell( er1, er2, er3, elvec )

! 可选：使用后清空指针，避免悬垂指针
nullify(er1, er2, er3)

内容的提问来源于stack exchange，提问作者Prokop Hapala