diff --git a/benchmarks/benchmark-dgemm.C b/benchmarks/benchmark-dgemm.C index 3867ed106..f2a37ecf5 100644 --- a/benchmarks/benchmark-dgemm.C +++ b/benchmarks/benchmark-dgemm.C @@ -40,9 +40,11 @@ typedef FFLAS::Timer TTimer; #ifndef __SGEMM__ typedef double Floats; #define CBLAS_GEMM cblas_dgemm +#define CUBLAS_GEMM cublasDgemm #else typedef float Floats; #define CBLAS_GEMM cblas_sgemm +#define CUBLAS_GEMM cublasDgemm #endif @@ -68,7 +70,7 @@ int main(int argc, char** argv) { FFLAS::parseArguments(argc,argv,as); - typedef Givaro::ModularBalanced Field; + typedef Givaro::Modular Field; typedef Field::Element Element; Field F(q); @@ -77,7 +79,7 @@ int main(int argc, char** argv) { double time=0.0;// time2=0.0; Element * A, * B, * C; - +/* if (iter>1) { if (!file1.empty()){ FFLAS::ReadMatrix (file1.c_str(),F,n,n,A); @@ -105,14 +107,19 @@ int main(int argc, char** argv) { C = FFLAS::fflas_new(n*n); + #if defined(CUDA_BLAS) + CUBLAS_GEMM ('n', 'n', n,n,n, F.one, + A, n, B, n, F.zero, C,n); // @fixme CUBLAS ALWAYS COLUMN MAJOR + #else CBLAS_GEMM (CblasRowMajor, CblasNoTrans, CblasNoTrans, n,n,n, F.one, A, n, B, n, F.zero, C,n); + #endif FFLAS::fflas_delete( A); FFLAS::fflas_delete( B); FFLAS::fflas_delete( C); } - +*/ for (size_t it=0;it(n*n); chrono.clear(); chrono.start(); + + #if defined(CUDA_BLAS) + + // Allocate device storage for A,B,C + Element *d_A, *d_B, *d_C; + cudaMalloc((void**)&d_A, n*n*sizeof(Element)); + cudaMalloc((void**)&d_B, n*n*sizeof(Element)); + cudaMalloc((void**)&d_C, n*n*sizeof(Element)); + + // Create cublas instance + cublasHandle_t handle; + cublasCreate(&handle); + + // Set input matrices on device + cublasSetMatrix(n, n, sizeof(Element), A, n, d_A, n); + cublasSetMatrix(n, n, sizeof(Element), B, n, d_B, n); + cublasSetMatrix(n, n, sizeof(Element), C, n, d_C, n); + + +// CUBLAS_GEMM ( 'n', 'n', n,n,n, F.one, A, n, B, n, F.zero, C,n); +CUBLAS_GEMM (handle, CUBLAS_OP_N, CUBLAS_OP_N, n,n,n, &F.one, d_A, n, d_B, n, &F.zero, d_C,n); + + // Retrieve result matrix from device + cublasGetMatrix(n, n, sizeof(Element), d_C, n, C, n); + + #else CBLAS_GEMM (CblasRowMajor, CblasNoTrans, CblasNoTrans, n,n,n, F.one, A, n, B, n, F.zero, C,n); + #endif + + std::cout << "C = A * B: " << fmod(C[0], q) << " expected " << fmod(A[0] * B[0], q) << std::endl; + chrono.stop(); time+=chrono.usertime(); diff --git a/configure.ac b/configure.ac index 04df5b84f..9f228dd49 100644 --- a/configure.ac +++ b/configure.ac @@ -287,7 +287,7 @@ FF_OPENBLAS_NUM_THREADS # BLAS_LIBS="-L/${BLAS_PATH} ${LAPACK_LIBS} ${BLAS_LIBS}" # AC_SUBST(BLAS_LIBS) -# FF_CHECK_CUDA + FF_CHECK_CUDA # AM_CONDITIONAL(FFLASFFPACK_HAVE_BLAS, test "x$BLAS_FOUND" != "xfalse") diff --git a/fflas-ffpack/config-blas.h b/fflas-ffpack/config-blas.h index a38023837..86aa91126 100644 --- a/fflas-ffpack/config-blas.h +++ b/fflas-ffpack/config-blas.h @@ -42,9 +42,14 @@ #ifdef __FFLASFFPACK_HAVE_MKL #include - #endif +#ifdef __FFLASFFPACK_HAVE_CUDA +#define __FFLASFFPACK_HAVE_CUBLAS +#include +#include +#include +#endif #ifndef CBLAS_INT #ifdef blasint /* openblas */ @@ -56,16 +61,6 @@ #endif /* blasint */ #endif /* CBLAS_INT */ -#ifdef CUDA_BLAS - -#define sgemv_ cublas_sgemv -#define sgemm_ cublas_sgemm -#define strsm_ cublas_strsm -#define strmm_ cublas_strmm - -#endif // CUDA_BLAS - - #ifndef __FFLASFFPACK_HAVE_MKL #define CBLAS_ENUM_DEFINED_H diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl index 334046af7..1f59e1337 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl @@ -41,6 +41,7 @@ #include "fflas-ffpack/fflas/fflas_igemm/igemm.h" #endif +#include "fflas-ffpack/utils/fflas_io.h" namespace FFLAS { // F is a field supporting delayed reductions @@ -249,12 +250,41 @@ namespace FFLAS { FFLASFFPACK_check(lda); FFLASFFPACK_check(ldb); FFLASFFPACK_check(ldc); -#if defined(__FFLASFFPACK_OPENBLAS_NUM_THREADS) and not defined (__FFLASFFPACK_OPENBLAS_NT_ALREADY_SET) - openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); -#endif - cblas_dgemm (CblasRowMajor, (CBLAS_TRANSPOSE) ta, (CBLAS_TRANSPOSE) tb, - (int)m, (int)n, (int)k, (Givaro::DoubleDomain::Element) alpha, - Ad, (int)lda, Bd, (int)ldb, (Givaro::DoubleDomain::Element) beta, Cd, (int)ldc); + + #if defined(__FFLASFFPACK_HAVE_CUBLAS) + // Allocate device storage for A,B,C + using Element = Givaro::DoubleDomain::Element; + Element *d_A, *d_B, *d_C; + cudaMalloc((void**)&d_A, m*k*sizeof(Element)); + cudaMalloc((void**)&d_B, k*n*sizeof(Element)); + cudaMalloc((void**)&d_C, m*n*sizeof(Element)); + + // Create cublas instance + cublasHandle_t handle; + cublasCreate(&handle); + + // Set input matrices on device + cublasSetMatrix(k, m, sizeof(Element), Ad, lda, d_A, k); // @note Device's leading dimensions are for column major matrices. + cublasSetMatrix(n, k, sizeof(Element), Bd, ldb, d_B, n); + cublasSetMatrix(n, m, sizeof(Element), Cd, ldc, d_C, n); + + // @note As cublas works in column-major, we inverse the transpose here. + auto d_ta = (ta == FFLAS_TRANSPOSE::FflasNoTrans) ? CUBLAS_OP_N : CUBLAS_OP_T; + auto d_tb = (tb == FFLAS_TRANSPOSE::FflasNoTrans) ? CUBLAS_OP_N : CUBLAS_OP_T; + + cublasDgemm(handle, d_tb, d_ta, n,m,k, &alpha, d_B, n, d_A, k, &beta, d_C, n); + + // Retrieve result matrix from device + cublasGetMatrix(n, m, sizeof(Element), d_C, n, Cd, ldc); + #else + #if defined(__FFLASFFPACK_OPENBLAS_NUM_THREADS) and not defined (__FFLASFFPACK_OPENBLAS_NT_ALREADY_SET) + openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); + #endif + + cblas_dgemm (CblasRowMajor, (CBLAS_TRANSPOSE) ta, (CBLAS_TRANSPOSE) tb, + (int)m, (int)n, (int)k, (Givaro::DoubleDomain::Element) alpha, + Ad, (int)lda, Bd, (int)ldb, (Givaro::DoubleDomain::Element) beta, Cd, (int)ldc); + #endif } inline void fgemm (const Givaro::FloatDomain& F, @@ -272,12 +302,41 @@ namespace FFLAS { FFLASFFPACK_check(ldb); FFLASFFPACK_check(ldc); -#if defined(__FFLASFFPACK_OPENBLAS_NUM_THREADS) and not defined (__FFLASFFPACK_OPENBLAS_NT_ALREADY_SET) - openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); -#endif - cblas_sgemm (CblasRowMajor, (CBLAS_TRANSPOSE) ta, (CBLAS_TRANSPOSE) tb, - (int)m, (int)n, (int)k, (Givaro::FloatDomain::Element) alpha, - Ad, (int)lda, Bd, (int)ldb, (Givaro::FloatDomain::Element) beta,Cd, (int)ldc); + #if defined(__FFLASFFPACK_HAVE_CUBLAS) + // Allocate device storage for A,B,C + using Element = Givaro::FloatDomain::Element; + Element *d_A, *d_B, *d_C; + cudaMalloc((void**)&d_A, m*k*sizeof(Element)); + cudaMalloc((void**)&d_B, k*n*sizeof(Element)); + cudaMalloc((void**)&d_C, m*n*sizeof(Element)); + + // Create cublas instance + cublasHandle_t handle; + cublasCreate(&handle); + + // Set input matrices on device + cublasSetMatrix(k, m, sizeof(Element), Ad, lda, d_A, k); // @note Device's leading dimensions are for column major matrices. + cublasSetMatrix(n, k, sizeof(Element), Bd, ldb, d_B, n); + cublasSetMatrix(n, m, sizeof(Element), Cd, ldc, d_C, n); + + // @note As cublas works in column-major, we inverse the transpose here. + auto d_ta = (ta == FFLAS_TRANSPOSE::FflasNoTrans) ? CUBLAS_OP_N : CUBLAS_OP_T; + auto d_tb = (tb == FFLAS_TRANSPOSE::FflasNoTrans) ? CUBLAS_OP_N : CUBLAS_OP_T; + + cublasSgemm(handle, d_tb, d_ta, n,m,k, &alpha, d_B, n, d_A, k, &beta, d_C, n); + + // Retrieve result matrix from device + cublasGetMatrix(n, m, sizeof(Element), d_C, n, Cd, ldc); + #else + #if defined(__FFLASFFPACK_OPENBLAS_NUM_THREADS) and not defined (__FFLASFFPACK_OPENBLAS_NT_ALREADY_SET) + openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); + #endif + cblas_sgemm (CblasRowMajor, (CBLAS_TRANSPOSE) ta, (CBLAS_TRANSPOSE) tb, + (int)m, (int)n, (int)k, (Givaro::FloatDomain::Element) alpha, + Ad, (int)lda, Bd, (int)ldb, (Givaro::FloatDomain::Element) beta,Cd, (int)ldc); + #endif + + } inline void fgemm (const Givaro::ZRing& F, diff --git a/macros/cuda-check.m4 b/macros/cuda-check.m4 index 7fc760932..1a391d087 100644 --- a/macros/cuda-check.m4 +++ b/macros/cuda-check.m4 @@ -63,54 +63,53 @@ do if test -r "$CUDA_HOME/include/cuda.h" ; then CUDA_CFLAGS="-I${CUDA_HOME}/include" CUDA_PATH="-L${CUDA_HOME}/lib64" - CUDA_LIBS="-L${CUDA_HOME}/lib64 -lcusparse" + CUDA_LIBS="-L${CUDA_HOME}/lib64 -lcusparse -lcublas -lcudart" else echo "($CUDA_HOME) seems an invalid CUDA prefix" echo "Searching CUDA in PATH" CUDA_CFLAGS="" - CUDA_LIBS="-lcusparse" + CUDA_LIBS="-lcusparse -lcublas -lcudart" fi else CUDA_CFLAGS="" - CUDA_LIBS="-lcusparse" + CUDA_LIBS="-lcusparse -lcublas -lcudart" fi CXXFLAGS="${CXXFLAGS} ${CUDA_CFLAGS}" LIBS="${LIBS} ${CUDA_LIBS}" CODE_CUDA=`cat macros/CodeChunk/cuda.C` + dnl This checks if cuda is installed, which ensures that cublas is available. + AC_TRY_LINK( [ - #include + #include ], - [ CUresult a;], + [ CUresult a; ], [ - dnl # See if we are running CUDA 4.0 with --enable-cxx - AC_TRY_RUN( - [ ${CODE_CUDA} ], - [ - AC_MSG_RESULT(found) - AC_DEFINE(HAVE_CUDA,1,[Define if CUDA is installed]) - - dnl CUDA_VERSION="" dnl I could find it but why is it here ? - CUDA_LIBS="${CUDA_PATH} -lcusparse" - dnl AC_SUBST(CUDA_VERSION) - AC_SUBST(CUDA_LIBS) - AC_SUBST(CUDA_CFLAGS) - break; - ],[ - AC_MSG_RESULT(no : cuda is too old or not found) - dnl AC_SUBST(CUDA_VERSION) - ],[ dnl This should never happen - AC_MSG_RESULT(no) - ]) + dnl # See if we are running CUDA 4.0 with --enable-cxx + AC_TRY_RUN( + [ ${CODE_CUDA} ], + [ + AC_MSG_RESULT(found) + + AC_SUBST(CUDA_LIBS) + AC_SUBST(CUDA_CFLAGS) + AC_DEFINE(HAVE_CUDA,1,[Define if CUDA blas is installed]) + break; + ],[ + AC_MSG_RESULT(no : cuda is too old or not found) + ],[ + AC_MSG_RESULT(no) + ] + ) ],[ AC_MSG_RESULT(unknown) echo "WARNING: You appear to be cross compiling, so there is no way to determine" echo "whether your CUDA version is new enough. I am assuming it is." AC_SUBST(CUDA_CFLAGS) AC_SUBST(CUDA_LIBS) - AC_DEFINE(HAVE_CUDA,1,[Define if CUDA is installed]) + AC_DEFINE(HAVE_CUDA,1,[Define if CUDA blas is installed]) ]) unset CUDA_CFLAGS unset CUDA_LIBS