在Fortran中传递一个不连续的数组部分

Question

我使用intel fortran编译器和intel mkl进行性能检查.我正在通过类似的调用将一些数组部分传递给Fortran 77接口

call dgemm( transa,transb,sz_s,P,P,&
            a, Ts_tilde,&
            sz_s,R_alpha,P,b,tr(:sz_s,:),sz_s)

很明显,tr(:sz_s,:)在内存中是不连续的,并且Fortran 77接口期望连续阻塞并为此创建一个临时块.

我想知道的是,如果我在tr的代码中显式创建临时数组并在操作之前和之后来回复制来自该临时数据的信息会有所不同,或者与创建临时的编译器本身相同从表现的角度来看？我想编译器总是会更有效率.

当然,我们欢迎任何有关消除这些临时工作的建议.

还有一点,如果我显然使用库的Fortran 95接口,对更简单的测试问题进行类似的调用,则不会发出创建临时的警告.然后我在mkl的手册中读到Fortran 95接口使用假定的形状数组,这解释了为什么不创建临时数据.

但是在那时,我似乎无法使用一些支持函数,如计时例程.也就是说,intel mkl有一些定时支持函数,但是如果我将它们与mkl_service例程一起使用,那么我得到'这个名字没有类型,并且dsecnd必须有一个显式类型'错误.对这个问题的任何想法也是受欢迎的.一个简单的例子给出了

program dgemm95_test
! some modules for Fortran 95 interface
use mkl_service
use mkl95_precision
use mkl95_blas
!
implicit none
!
double precision, dimension(4,3) :: a
double precision, dimension(6,4) :: b
double precision, dimension(5,5) :: r ! result array
double precision, dimension(3,2) :: dummy_b
!
character(len=1) :: transa
character(len=1) :: transb
!
double precision :: alpha, beta, t1, t2, t
integer :: sz1, sz2

! initialize some variables
alpha = 1.0
beta = 0.0
a = 2.3
b = 4.5
r = 0.0
transa = 'n'
transb = 'n'
dummy_b = 0.0
! Fortran 95 interface
t1 = dsecnd()
call gemm( a, b(4:6,1:3:2), r(2:5,3:4),&
 transa, transb, alpha, beta )
t2 = dsecnd()
!
write(*,*) r
dummy_b  = r(2:4,4:5)
!
end program dgemm95_test

Answer 1

将数组部分传递给假定大小的数组伪参数(旧例程使用)时,临时是绝对必要的,因为数组部分在内存中不连续.

你当然可以制作自己的临时数组.它是否会更快取决于许多因素.其中重要的是临时是在堆栈上还是在堆上分配.英特尔Fortran编译器兼具两者,有编译器开关来控制行为(-heap-arrays n),它可以依赖于数组大小.堆栈分配要快得多,通常是默认值.默认情况下,您可能用于自己的临时的自动数组也会在堆栈上分配.小心堆栈上的大型数组,您可以轻松地溢出它并导致崩溃.

我建议你做一个性能测试,如果它不是太慢,可以使用更简单的变体.可能它将是Fortran 95界面,但你应该测量时间,真的.

至于时间,MKL手册页second()//dsecnd()你必须包括mkl_lapack.fi但不谈论任何Fortran95接口.你可以逃避声明外部双精度,但我会使用包含.或者system_clock()用作便携式标准Fortran 95.