【发布时间】：2014-04-11 00:59:49
【问题描述】：
在英特尔 Fortran 等现代 Fortran 编译器中是否可以在运行时确定数组步长？例如，我可能想对数组部分执行快速傅里叶变换 (FFT)：
program main
complex(8),allocatable::array(:,:)
allocate(array(17, 17))
array = 1.0d0
call fft(array(1:16,1:16))
contains
subroutine fft(a)
use mkl_dfti
implicit none
complex(8),intent(inout)::a(:,:)
type(dfti_descriptor),pointer::desc
integer::stat
stat = DftiCreateDescriptor(desc, DFTI_DOUBLE, DFTI_COMPLEX, 2, shape(a) )
stat = DftiCommitDescriptor(desc)
stat = DftiComputeForward(desc, a(:,1))
stat = DftiFreeDescriptor(desc)
end subroutine
end program
但是，需要明确告知 MKL Dfti* 例程的数组步幅。 翻阅参考手册，我没有发现任何返回步幅信息的内在函数。 一些有趣的资源是here 和here，它们讨论了是否复制数组部分以及英特尔 Fortran 如何在内部处理数组。 我宁愿不局限于英特尔目前使用其数组描述符的方式。
如何计算步幅信息？请注意，通常我希望fft 例程（或任何类似例程）不需要传入有关数组的任何其他信息。
编辑：
我已经验证了在这种情况下不是创建了一个临时数组，这是我在 Intel(R) Visual Fortran Compiler XE 14.0.2.176 [Intel( R) 64]，禁用优化并将堆数组设置为 0。
program main
implicit none
real(8),allocatable::a(:,:)
pause
allocate(a(8192,8192))
pause
call random_number(a)
pause
call foo(a(:4096,:4096))
pause
contains
subroutine foo(a)
implicit none
real(8)::a(:,:)
open(unit=16, file='a_sum.txt')
write(16, *) sum(a)
close(16)
end subroutine
end program
监控内存使用情况，很明显永远不会创建一个临时数组。
编辑 2：
module m_foo
implicit none
contains
subroutine foo(a)
implicit none
real(8),contiguous::a(:,:)
integer::i, j
open(unit=16, file='a_sum.txt')
write(16, *) sum(a)
close(16)
call nointerface(a)
end subroutine
end module
subroutine nointerface(a)
implicit none
real(8)::a(*)
end subroutine
program main
use m_foo
implicit none
integer,parameter::N = 8192
real(8),allocatable::a(:,:)
integer::i, j
real(8)::count
pause
allocate(a(N, N))
pause
call random_number(a)
pause
call foo(a(:N/2,:N/2))
pause
end program
编辑 3：
该示例说明了我想要实现的目标。有一个 16x16 的连续数组，但我只想转换上面的 4x4 数组。第一次调用只是传入数组部分，但它不会在数组的左上角返回一个。第二个调用设置适当的步幅，a 随后包含正确的上 4x4 数组。上层 4x4 数组相对于整个 16x16 数组的步幅不是一个。
program main
implicit none
complex(8),allocatable::a(:,:)
allocate(a(16,16))
a = 0.0d0
a(1:4,1:4) = 1.0d0
call fft(a(1:4,1:4))
write(*,*) a(1:4,1:4)
pause
a = 0.0d0
a(1:4,1:4) = 1.0d0
call fft_stride(a(1:4,1:4), 1, 16)
write(*,*) a(1:4,1:4)
pause
contains
subroutine fft(a) !{{{
use mkl_dfti
implicit none
complex(8),intent(inout)::a(:,:)
type(dfti_descriptor),pointer::desc
integer::stat
stat = DftiCreateDescriptor(desc, DFTI_DOUBLE, DFTI_COMPLEX, 2, shape(a) )
stat = DftiCommitDescriptor(desc)
stat = DftiComputeForward(desc, a(:,1))
stat = DftiFreeDescriptor(desc)
end subroutine !}}}
subroutine fft_stride(a, s1, s2) !{{{
use mkl_dfti
implicit none
complex(8),intent(inout)::a(:,:)
integer::s1, s2
type(dfti_descriptor),pointer::desc
integer::stat
integer::strides(3)
strides = [0, s1, s2]
stat = DftiCreateDescriptor(desc, DFTI_DOUBLE, DFTI_COMPLEX, 2, shape(a) )
stat = DftiSetValue(desc, DFTI_INPUT_STRIDES, strides)
stat = DftiCommitDescriptor(desc)
stat = DftiComputeForward(desc, a(:,1))
stat = DftiFreeDescriptor(desc)
end subroutine !}}}
end program
【问题讨论】：
-
我怀疑，当 MKL FFT 例程谈论步幅时，他们真正关心的是您正在执行的 FFT 的数据类型和类型，即您是否正在对真实或复杂数据执行 ND FFT。
-
根据文档，
DftiCreateDescriptor将使用无填充的默认假设，如果您最终将假设形状传递给低级例程，则会发生这种情况。如果你想传递一块array，那么传递整个数组并设置步幅将比传递数组切片更有效，如果它不连续，则会创建一个临时数组。 -
您确定它会临时创建一个数组吗？我想如果有一个明确的接口就没有临时的：nf.nci.org.au/facilities/software/FORTRAN/Intel10/doc/main_for/… 见右下角。
-
@Yossarian 实际上他们确实允许输入和输出内存以非常通用的方式布局，只要有规律的步幅。我已经使用 C 接口确认了这一点。
-
@bdforbes 一旦您进入 mkl 例程，就会创建临时数组，当然不在您的代码中，因为您的虚拟参数是一个假定的形状。