SUBROUTINE SGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
      
      USE ISO_C_BINDING
      IMPLICIT NONE
      
!     .. Scalar Arguments ..
      REAL ALPHA,BETA
      INTEGER K,LDA,LDB,LDC,M,N
      CHARACTER(LEN=*) TRANSA,TRANSB
      CHARACTER(KIND=c_char) CTA, CTB
!     ..
!     .. Array Arguments ..
      REAL A(LDA,*),B(LDB,*),C(LDC,*)
! Define the INTERFACE to the NVIDIA C code cublasSgemm.
! This version of SGEMM is used in a user application
! that calls LAPACK single precision routines, or makes
! other uses of that code.
 
      INTERFACE
! This is what the NVIDIA code expects for its inputs:      
!      void cublasSgemm (char transa, char transb, int m, int n,
!         int k, float alpha, const float *A, int lda,
!         const float *B, int ldb, float beta,
!         float *C, int ldc)
          subroutine c_sgemm(transa, transb, m, n, k,&
             alpha, A, lda, B, ldb, beta, c, ldc)bind(C,name='cublasSgemm')
             USE, INTRINSIC :: iso_c_binding, ONLY : c_int, c_float, c_char 
             
             character(KIND=c_char),value :: transa, transb 
             integer(c_int),value :: m,n,k,lda,ldb,ldc
             real(c_float),value :: alpha,beta
             real(c_float) :: A(lda,*),B(ldb,*),C(ldc,*)
          end subroutine c_sgemm
      END INTERFACE 

! The calculation, excepting initialization and finalization,
! is done with the NVIDIA C routine 'cublasSgemm.'
! A local name c_sgemm is used in Fortran.
! The name c_sgemm could be replaced by NVIDIA's name if one chose. 
      cta = transa(1:1)
      ctb = transb(1:1)
      call c_sgemm(cta, ctb, &
       m, n, k, alpha, A, lda, B, ldb, beta, c, ldc)
      return     
!     ..
      END SUBROUTINE SGEMM