在MPI程序中调用cublas和cula是很经常的事情，然而很诡异的是，利用多块GPU卡做MPI并行计算，我却发现，多块卡并行时，比单块卡计算的时间要长得多。为了查清楚原因所在，写了如下简单的程序来做单GPU计算和双GPU计算(调用zgesv，用来做矩阵求逆）：

PROGRAM cula_test

use cudafor

use cula_status

use cula_lapack

use cula_lapack_device_pgfortran

IMPLICIT NONE

include 'mpif.h'

INTEGER :: n

complex*16, allocatable::A(:,:),U(:,:)

integer,allocatable::ipiv(:)

integer I,J,info,MPIerror,node,Nnodes

real*8 c,d

real*4 t1,t2

external cula_initialize

external cula_shutdown

external cudasetdevice

call MPI_Init( MPIerror )

call MPI_Comm_Rank( MPI_Comm_World, Node, MPIerror )

call MPI_Comm_Size( MPI_Comm_World, Nnodes, MPIerror )

if(node.eq.0) info=cudasetdevice(0)

if(node.eq.1) info=cudasetdevice(1)

info = cula_initialize()

n = 3000

ALLOCATE(A(n,n),U(n,n), ipiv(n))

do I = 1,N

  do J = 1, N

     call random_number(c)

     call random_number(d)

     A(I,J)=dcmplx(c,d)

  enddo

enddo

U(:,:)=(0.d0,0.d0)

do I = 1, N

  U(I,I)=(1.d0,0.d0)

enddo

call cpu_time(t1)

info= cula_zgesv(n,n,A,n,ipiv,U,n)

call cpu_time(t2)

print *,'GPU: ', U(1,1),t2-t1

deallocate(A,U,ipiv)

call cula_shutdown()

call MPI_FINALIZE(MPIerror)

END

当用一个进程，一块K20c卡时，输出结果是：

GPU:   (-8.6050432536450713E-002,-0.1393513431401034)   0.7435620

当用两个进程，两块K20c卡时，输出结果是：

GPU:   (-8.6050432536450713E-002,-0.1393513431401034)    7.932089

GPU:   (-8.6050432536450713E-002,-0.1393513431401034)    8.396360

很明显的，当用两个进程时，我只是让两块卡做了完全相同的事情（对3000×3000的矩阵求逆），调用了同一个cula的子程序，两个进程之间不存在任何直接的数据通讯，彼此完全独立，但所需的时间却是只用一块卡时的十倍。我也曾经做测试过2000×2000的矩阵，所需的时间相差20倍。实在感觉匪夷所思。

接下来我测试了cublas中的zgemm程序。源代码如下：

PROGRAM cublas_test

use cudafor

use cublas_pgfortran

IMPLICIT NONE

include 'mpif.h'

INTEGER :: n

complex*16, allocatable::A(:,:),B(:,:),C(:,:)

integer info, I,J,MPIerror,node,Nnodes

real*8 e,f

real*4 t1,t2

character*1 TRANSA, TRANSB

complex*16 ALPHA, BETA

external cudasetdevice

ALPHA=(1.D0,0.D0)

BETA=(0.D0,0.D0)

TRANSA='N'

TRANSB='N'

call MPI_Init( MPIerror )

call MPI_Comm_Rank( MPI_Comm_World, Node, MPIerror )

call MPI_Comm_Size( MPI_Comm_World, Nnodes, MPIerror )

if(node.eq.0) info=cudasetdevice(0)

if(node.eq.1) info=cudasetdevice(1)

N = 8000  !matrix size

ALLOCATE(A(n,n),B(n,n), C(n,n))

do I = 1,N

  do J = 1, N

     call random_number(e)

     call random_number(f)

     A(I,J)=dcmplx(e,f)

  enddo

enddo

B = A

C = A

call cpu_time(t1)

call  cublas_zgemm(TRANSA,TRANSB,n,n,n,ALPHA,A,n,n,B,n,n,BETA,C,n,n)

call cpu_time(t2)

print *,'GPU: ', C(1,1),t2-t1

deallocate(A,B,C)

call MPI_FINALIZE(MPIerror)

END

当用一个进程，一块K20c卡时，输出结果是：

GPU:   (-2.032549874426836,4008.562512905637)    5.963408

当用两个进程，两块K20c卡时，输出结果是：

GPU:   (-2.032549874426836,4008.562512905637)    5.928651

GPU:   (-2.032549874426836,4008.562512905637)    5.944520

我们看到，调用cublas中的函数时，不论是用一块卡，还是并行两块卡时，用的时间基本是一样的。

我只是测试了zgesv（cula）和zgemm（cublas）两个子程序，得到的结论是，调用cublas时，多块显卡并行时是没有问题的，而调用cula时，多块卡并行（不管进程间有没有数据通讯）时会使得时间变得不可思议的长。我也曾经试过，让一块卡运行完毕，再让另一块开始运行，则各自所花时间又变得正常。