#include <stdlib.h> 

#include <stdio.h> 

#include <time.h> 

#include <sys/time.h> 

#include "cublas.h" 

#include "cuda.h" 

//#define N_SIZES 8 

#define N_SIZES 1 

//int sizes[] = { 256, 512, 768, 1024, 1280, 1536, 1792, 2048 }; 

int sizes[] = { 1536 }; 

//#define N_TRIALS 500 

#define N_TRIALS 1 

//#define NI_TRIALS 5000 

#define NI_TRIALS 5000 

void print_dmatrix(double *matrix, int nrows, int ncolumns) { 

  int m, n; 

  // print result 

  for (m = 0; m < nrows; m++) { 

    for (n = 0; n < ncolumns; n++) { 

      printf("%f ", matrix[m * ncolumns + n]); 

    } 

    printf("\n"); 

  } 

} 

int main(int argc, void *argv[]) { 

  float           *hA_s, *hb_s, *dA_s, *db_s, *dresult_s, *hresult_s; 

  double          *hA_d, *hb_d, *dA_d, *db_d, *dresult_d, *hresult_d; 

  cuComplex       *hA_c, *hb_c, *dA_c, *db_c, *dresult_c, *hresult_c; 

  cuDoubleComplex *hA_z, *hb_z, *dA_z, *db_z, *dresult_z, *hresult_z; 

  int M, N, K;  // sizes of matrices 

  int s, i, t;  // counters 

  cublasStatus status; 

  cudaError_t error; 

  timeval start_time, end_time, elapsed_time; 

  double single_times[N_SIZES]; 

  double double_times[N_SIZES]; 

  double complex_times[N_SIZES]; 

  double dcomplex_times[N_SIZES]; 

  srand48(time(NULL)); 

  // initialize CUBLAS 

  status = cublasInit(); 

  if (status != CUBLAS_STATUS_SUCCESS) printf("Error initializing CUBLAS!\n"); 

  for (s = 0; s < N_SIZES; s++) { 

    single_times[s] = 0.0; 

    double_times[s] = 0.0; 

    complex_times[s] = 0.0; 

    dcomplex_times[s] = 0.0; 

  } 

  printf("Running timing for cuBLAS\n"); 

  for (t = 0; t < N_TRIALS; t++ ) { 

    //    if ((t + 1) % 10 == 0) { 

    //      printf("."); 

    //      fflush(NULL); 

    //    } 

    for (s = 0; s < N_SIZES; s++) { 

      M = sizes[s]; 

      N = sizes[s]; 

      K = sizes[s]; 

      // start float timing 

      printf("Starting float timing...\n"); 

      // allocate host memory 

      hA_s = (float *) malloc(M * K * sizeof(float)); 

      hb_s = (float *) malloc(K * N * sizeof(float)); 

      hresult_s = (float *) malloc(M * N * sizeof(float)); 

      // allocate device memory 

      error = cudaMalloc(&dA_s, M * K * sizeof(float)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&db_s, K * N * sizeof(float)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&dresult_s, M * N * sizeof(float)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      // initialize host memory 

      for (i = 0; i < M * K; i++) hA_s[i] = drand48(); 

      for (i = 0; i < K * N; i++) hb_s[i] = drand48(); 

      // transfer to GPU 

      error = cudaMemcpy(dA_s, hA_s, M * K * sizeof(float), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      error = cudaMemcpy(db_s, hb_s, K * N * sizeof(float), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      // perform multiplication 

      /*  

	 C = alpha * op(A) * op(B) + beta * C 

	 where op(X) = X or op(X) = transpose(X) 

	 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)   

      */ 

      gettimeofday(&start_time, NULL); 

      for (t = 0; t < NI_TRIALS; t++ ) { 

	cublasSgemm('N', 'N',  

		    M, N, K,   

		    1.0, dA_s, M, 

		    db_s, K, 

		    0.0, dresult_s, M); 

      } 

      gettimeofday(&end_time, NULL); 

      timersub(&end_time, &start_time, &elapsed_time); 

      single_times[s] += elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0; 

      // transfer result to CPU 

      error = cudaMemcpy(hresult_s, dresult_s, M * N * sizeof(float), cudaMemcpyDeviceToHost); 

      if (error != cudaSuccess) printf("Error copying memory, device to host!\n"); 

      /* 

	print_dmatrix(hsA, M, K); 

	printf("--\n"); 

	print_dmatrix(hsb, K, N); 

	printf("--\n"); 

	print_dmatrix(hsresult, M, N); 

	printf("\n"); 

      */ 

      free(hA_s); 

      free(hb_s); 

      free(hresult_s); 

      cudaFree(dA_s); 

      cudaFree(db_s); 

      cudaFree(dresult_s); 

      // start double timing 

      printf("Starting double timing...\n"); 

      // allocate host memory 

      hA_d = (double *) malloc(M * K * sizeof(double)); 

      hb_d = (double *) malloc(K * N * sizeof(double)); 

      hresult_d = (double *) malloc(M * N * sizeof(double)); 

      // allocate device memory 

      error = cudaMalloc(&dA_d, M * K * sizeof(double)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&db_d, K * N * sizeof(double)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&dresult_d, M * N * sizeof(double)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      // initialize host memory 

      for (i = 0; i < M * K; i++) hA_d[i] = drand48(); 

      for (i = 0; i < K * N; i++) hb_d[i] = drand48(); 

      // transfer to GPU 

      error = cudaMemcpy(dA_d, hA_d, M * K * sizeof(double), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      error = cudaMemcpy(db_d, hb_d, K * N * sizeof(double), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      // perform multiplication 

      /*  

	 C = alpha * op(A) * op(B) + beta * C 

	 where op(X) = X or op(X) = transpose(X) 

	 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)   

      */ 

      gettimeofday(&start_time, NULL); 

      for (t = 0; t < NI_TRIALS; t++ ) { 

	cublasDgemm('N', 'N',  

		    M, N, K,   

		    1.0, dA_d, M, 

		    db_d, K, 

		    0.0, dresult_d, M); 

      } 

      gettimeofday(&end_time, NULL); 

      timersub(&end_time, &start_time, &elapsed_time); 

      //printf("time = %f seconds\n", elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0); 

      double_times[s] += elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0; 

      // transfer result to CPU 

      error = cudaMemcpy(hresult_d, dresult_d, M * N * sizeof(double), cudaMemcpyDeviceToHost); 

      if (error != cudaSuccess) printf("Error copying memory, device to host!\n"); 

      /* 

	print_dmatrix(hA, M, K); 

	printf("--\n"); 

	print_dmatrix(hb, K, N); 

	printf("--\n"); 

	print_dmatrix(hresult, M, N); 

	printf("\n"); 

      */ 

      free(hA_d); 

      free(hb_d); 

      free(hresult_d); 

      cudaFree(dA_d); 

      cudaFree(db_d); 

      cudaFree(dresult_d); 

      // start complex timing 

      printf("Starting complex timing...\n"); 

      // allocate host memory 

      hA_c = (cuComplex *) malloc(M * K * sizeof(cuComplex)); 

      hb_c = (cuComplex *) malloc(K * N * sizeof(cuComplex)); 

      hresult_c = (cuComplex *) malloc(M * N * sizeof(cuComplex)); 

      // allocate device memory 

      error = cudaMalloc(&dA_c, M * K * sizeof(cuComplex)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&db_c, K * N * sizeof(cuComplex)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&dresult_c, M * N * sizeof(cuComplex)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      // initialize host memory 

      for (i = 0; i < M * K; i++) {  

	hA_c[i].x = drand48(); 

	hA_c[i].y = drand48(); 

      } 

      for (i = 0; i < K * N; i++) { 

	hb_c[i].x = drand48(); 

	hb_c[i].y = drand48(); 

      } 

      // transfer to GPU 

      error = cudaMemcpy(dA_c, hA_c, M * K * sizeof(cuComplex), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      error = cudaMemcpy(db_c, hb_c, K * N * sizeof(cuComplex), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      // perform multiplication 

      /*  

	 C = alpha * op(A) * op(B) + beta * C 

	 where op(X) = X or op(X) = transpose(X) 

	 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)   

      */ 

      gettimeofday(&start_time, NULL); 

      for (t = 0; t < NI_TRIALS; t++ ) { 

	cublasCgemm('N', 'N',  

		    M, N, K,   

		    make_cuComplex(1.0, 0.0), dA_c, M, 

		    db_c, K, 

		    make_cuComplex(0.0, 0.0), dresult_c, M); 

      } 

      gettimeofday(&end_time, NULL); 

      timersub(&end_time, &start_time, &elapsed_time); 

      //printf("time = %f seconds\n", elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0); 

      complex_times[s] += elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0; 

      // transfer result to CPU 

      error = cudaMemcpy(hresult_c, dresult_c, M * N * sizeof(cuComplex), cudaMemcpyDeviceToHost); 

      if (error != cudaSuccess) printf("Error copying memory, device to host!\n"); 

      /* 

	print_dmatrix(hA, M, K); 

	printf("--\n"); 

	print_dmatrix(hb, K, N); 

	printf("--\n"); 

	print_dmatrix(hresult, M, N); 

	printf("\n"); 

      */ 

      free(hA_c); 

      free(hb_c); 

      free(hresult_c); 

      cudaFree(dA_c); 

      cudaFree(db_c); 

      cudaFree(dresult_c); 

      // start double complex timing 

      printf("Starting double complex timing...\n"); 

      // allocate host memory 

      hA_z = (cuDoubleComplex *) malloc(M * K * sizeof(cuDoubleComplex)); 

      hb_z = (cuDoubleComplex *) malloc(K * N * sizeof(cuDoubleComplex)); 

      hresult_z = (cuDoubleComplex *) malloc(M * N * sizeof(cuDoubleComplex)); 

      // allocate device memory 

      error = cudaMalloc(&dA_z, M * K * sizeof(cuDoubleComplex)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&db_z, K * N * sizeof(cuDoubleComplex)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      error = cudaMalloc(&dresult_z, M * N * sizeof(cuDoubleComplex)); 

      if (error != cudaSuccess) printf("Error allocating memory!\n"); 

      // initialize host memory 

      for (i = 0; i < M * K; i++) { 

	hA_z[i].x = drand48(); 

	hA_z[i].y = drand48(); 

      } 

      for (i = 0; i < K * N; i++) { 

	hb_z[i].x = drand48(); 

	hb_z[i].y = drand48(); 

      } 

      // transfer to GPU 

      error = cudaMemcpy(dA_z, hA_z, M * K * sizeof(cuDoubleComplex), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      error = cudaMemcpy(db_z, hb_z, K * N * sizeof(cuDoubleComplex), cudaMemcpyHostToDevice); 

      if (error != cudaSuccess) printf("Error copying memory, host to device!\n"); 

      // perform multiplication 

      /*  

	 C = alpha * op(A) * op(B) + beta * C 

	 where op(X) = X or op(X) = transpose(X) 

	 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)   

      */ 

      gettimeofday(&start_time, NULL); 

      for (t = 0; t < NI_TRIALS; t++ ) { 

	cublasZgemm('N', 'N',  

		    M, N, K,   

		    make_cuDoubleComplex(1.0, 0.0), dA_z, M, 

		    db_z, K, 

		    make_cuDoubleComplex(0.0, 0.0), dresult_z, M); 

      } 

      gettimeofday(&end_time, NULL); 

      timersub(&end_time, &start_time, &elapsed_time); 

      //printf("time = %f seconds\n", elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0); 

      dcomplex_times[s] += elapsed_time.tv_sec + elapsed_time.tv_usec / 1000000.0; 

      // transfer result to CPU 

      error = cudaMemcpy(hresult_z, dresult_z, M * N * sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost); 

      if (error != cudaSuccess) printf("Error copying memory, device to host!\n"); 

      /* 

	print_dmatrix(hA, M, K); 

	printf("--\n"); 

	print_dmatrix(hb, K, N); 

	printf("--\n"); 

	print_dmatrix(hresult, M, N); 

	printf("\n"); 

      */ 

      free(hA_z); 

      free(hb_z); 

      free(hresult_z); 

      cudaFree(dA_z); 

      cudaFree(db_z); 

      cudaFree(dresult_z); 

    } 

  } 

  // free resources 

  status = cublasShutdown(); 

  if (status != CUBLAS_STATUS_SUCCESS) printf("Error shutting down CUBLAS!\n"); 

  for (s = 0; s < N_SIZES; s++) { 

    double_times[s] = double_times[s] / (N_TRIALS * NI_TRIALS); 

    single_times[s] = single_times[s] / (N_TRIALS * NI_TRIALS); 

    complex_times[s] = complex_times[s] / (N_TRIALS * NI_TRIALS); 

    dcomplex_times[s] = dcomplex_times[s] / (N_TRIALS * NI_TRIALS); 

  } 

  printf("\nTrials: %d\n\n", N_TRIALS * NI_TRIALS); 

  for (s = 0; s < N_SIZES; s++) { 

    printf("Average time for multiplying %4d x %4d single matrices: %f sec\n",  

	   sizes[s], sizes[s], single_times[s]); 

    printf("GFLOPS: %f\n",  

	   2. * (sizes[s] / 1000.) * (sizes[s] / 1000.) * (sizes[s] / 1000.) / single_times[s]); 

  } 

  printf("\n"); 

  for (s = 0; s < N_SIZES; s++) { 

    printf("Average time for multiplying %4d x %4d double matrices: %f sec\n",  

	   sizes[s], sizes[s], double_times[s]); 

    printf("GFLOPS: %f\n",  

	   2. * (sizes[s] / 1000.) * (sizes[s] / 1000.) * (sizes[s] / 1000.) / double_times[s]); 

  } 

  for (s = 0; s < N_SIZES; s++) { 

    printf("Average time for multiplying %4d x %4d complex matrices: %f sec\n",  

	   sizes[s], sizes[s], complex_times[s]); 

    printf("GFLOPS: %f\n",  

	   8. * (sizes[s] / 1000.) * (sizes[s] / 1000.) * (sizes[s] / 1000.) / complex_times[s]); 

  } 

  for (s = 0; s < N_SIZES; s++) { 

    printf("Average time for multiplying %4d x %4d double complex matrices: %f sec\n",  

	   sizes[s], sizes[s], dcomplex_times[s]); 

    printf("GFLOPS: %f\n",  

	   8. * (sizes[s] / 1000.) * (sizes[s] / 1000.) * (sizes[s] / 1000.) / dcomplex_times[s]); 

  } 

}