PGI CUDA Fortran による CUBLAS 4.0 以降の使用方法

!
!      Copyright 2011, SofTek Systems, Inc.
!      All rights reserved.
!
program Sgemm_cublas

  use cudafor
  use cublas   ! for PGI 11.8 and later

  ! GPU device arrays
  real(4), device, allocatable, dimension(:,:) :: dA, dB, dC
  ! device properties
  type(cudaDeviceProp) :: prop
  ! device number
  integer idev
  ! command argument
  character*20 arg
  
  ! handle for CUBLAS
  type(cublasHandle) :: h1
  integer(4) :: iface = 0

  real(4), allocatable, dimension(:,:) :: a, b, c
  real(4) :: alpha = 1.0e0, beta  = 0.0e0
  ! for timing
  real(8) :: GFLOPS
  real(8) :: time
  integer(4) :: hz, clock0, clock1
  integer :: i, j, k

! BLAS interface
  print *, "Select BLAS interface? 0=traditional, 1=Portable legacy CUBLAS <4.0, 2=CUBLAS Version2"
  read (5,*) iface
  print *, "BLAS interface =", iface
  if (iface > 3) stop

! Matrix Size (N x N x N)
  print *, "What size of matrix : ?"
  read(5,*) n
  print *, "MatSize = ", n, " x ", n, " x ",n

! command argument and GPU device set
  nargs = command_argument_count()
  idev = 0
  do i = 1, nargs
    call get_command_argument(i,arg)
    if ((arg(1:7) .eq. "-device") .and. (i.lt.nargs)) then
      call get_command_argument(i+1,arg)
      read(arg,'(i2)') idev
    end if
  end do
  ! set a GPU device and print the properties
  istat = cudaSetDevice(idev)
  istat = cudaGetDeviceProperties(prop,idev)
  ilen = verify(prop%name, ' ', .true.)
  write (*,900) prop%name(1:ilen), &
                real(prop%clockRate)/1000.0, &
                real(prop%totalGlobalMem)/1024.0/1024.0

  ! arrays for HOST
  allocate(a(n,n), b(n,n), c(n,n))
  ! arrays for GPU
  allocate(dA(n,n), dB(n,n), dC(n,n))

  a = 1.0e0
  b = 5.5e0

  ! Initialize time function
  call system_clock(count_rate=hz)
  print *, "system_clock resolution: ", real(1.d0/hz)
  ! Start a clock
  call system_clock(count=clock0)

  ! for GPU
  dA = a
  dB = b
  if (beta .ne. 0.0d0) then
    dC = c
  endif

  if (iface == 0) then
    ! traditional BLAS interface
    call sgemm('n', 'n', n, n, n, alpha, dA, n, dB, n, beta, dC, n)

  else if (iface == 1) then
    ! Portable legacy CUBLAS interface
    call cublasSgemm('n', 'n', n, n, n, alpha, dA, n, dB, n, beta, dC, n)

  else if (iface ==2) then
    ! New CUBLAS 4.0 interface
    ! cublasOperation_t is defined as a integer namuber.
    ! "n" : CUBLAS_OP_N=0, "t" : CUBLAS_OP_T=1,  "c" : CUBLAS_OP_C=2
    h1 = cublasGetHandle()
    istat = cublasSgemm_v2(h1, 0, 0, n, n, n, alpha, dA, n, dB, n, beta, dC, n)

  end if

  ierr = cudathreadsynchronize()
  c = dC

  ! Stop a clock
  call system_clock(count=clock1)
  time= real((clock1-clock0)) / real (hz)

print *, "*** Verify results ***"

do j = 1, n
  do i = 1, n
    if (c(i,j) .ne. (5.50*real(n))) then
      print *, "error:  ",i,j,c(i,j)
    endif
  enddo
enddo

  GFLOPS = (real(n) * real(n) * real(n) * 2.0)
  print *, "Total Time (sec): ",time
  print *, "Total SGEMM GFLOPS: ",GFLOPS/(time*1.0d9)
  print *, ""

900 format('\nDevice:',a,', ',f6.1,' MHz clock, ',f6.1,' MB memory.\n')
end

●コンパイル・リンク メッセージ -Minfo 
[kato@photon29]$ pgfortran -o a.out Sgemm_cublas.cuf -fastsse -Mcuda -Minfo -lcublas
sgemm_cublas:
     48, Loop not vectorized/parallelized: contains call
     68, Memory set idiom, array assignment replaced by call to pgf90_mset4
     69, Memory set idiom, array assignment replaced by call to pgf90_mset4
    110, Loop not vectorized/parallelized: contains call