Advertisement
Guest User

Untitled

a guest
Dec 14th, 2012
64
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <cuda.h>
  2. #include <cuda_runtime_api.h>
  3. #include <stdio.h>
  4. #include <cublas_v2.h>
  5.  
  6. //const int R=16; // количество частотных интервалов всего
  7. //const float pi=3.1415926535897;
  8. //const float sigma=pi/R;
  9.  
  10. void mMul(int m, int n, int k, const float * a, const float * b, float * c)
  11. {
  12.     for (int i=0;i<m;i++)
  13.         for (int j = 0; j < k; ++j)
  14.         {
  15.             c[i*k+j]=0;
  16.             for (int l=0;l<n;l++)
  17.                 c[i*k+j]+=a[i*n+l]*b[l*k+j];
  18.         }
  19. }
  20.  
  21.  
  22. int main (int argc, char * argv[])
  23. {
  24.      //matrices A (m x n), B(n x k), C(m x k), D(k x l) and E (m x l)
  25.      const int m=11;
  26.      const int n=11;
  27.      const int k=11;
  28.      const int l=11;
  29.      
  30.      cudaError_t cudaStat;    
  31.      cublasStatus_t stat;
  32.      cublasHandle_t handle;
  33.      //initialization
  34.      
  35.      stat = cublasCreate(&handle);
  36.      if (stat != CUBLAS_STATUS_SUCCESS)
  37.      {
  38.       printf ("CUBLAS initialization failed\n");
  39.       return EXIT_FAILURE;
  40.      }
  41.      
  42.      //GPU pointers
  43.      float *dev_A, *dev_B, *dev_C, *dev_D, *dev_E;
  44.      //CPU pointers
  45.      float *host_A, *host_B, *host_C, *host_D, *host_E, *host_result;
  46.      
  47.      
  48.      //CPU memory allocation
  49.      //first matrix
  50.      host_A=new float [m*n];
  51.      
  52.      for (int i=0; i<m; i++)
  53.       for (int j=0; j<n; j++)
  54.            if (i==j)
  55.             host_A[i*n+j]=1;
  56.            else
  57.             host_A[i*n+j]=0;
  58.      
  59.      //the second
  60.      host_B=new float [n*k];
  61.      for (int i=0; i<n; i++)
  62.       for (int j=0; j<k; j++)
  63.            host_B[i*k+j]=i*k+j;
  64.      
  65.      //the third
  66.      host_D=new float [k*l];
  67.      for (int i=0; i<k; i++)
  68.       for (int j=0; j<l; j++)
  69.            host_D[i*l+j]=i*l+j;
  70.      
  71.      //result
  72.      host_C=new float [m*k];
  73.      host_E=new float [m*l];
  74.      
  75.      //GPU memory allocation
  76.      cudaStat = cudaMalloc ((void**)&dev_A, m*n*sizeof(float));
  77.      if (cudaStat != cudaSuccess) {
  78.       printf ("device memory allocation failed\n");
  79.       return EXIT_FAILURE;
  80.      }
  81.      cudaStat = cudaMalloc ((void**)&dev_B, n*k*sizeof(float));
  82.      if (cudaStat != cudaSuccess) {
  83.       printf ("device memory allocation failed\n");
  84.       return EXIT_FAILURE;
  85.      }
  86.      cudaStat = cudaMalloc ((void**)&dev_C, m*k*sizeof(float));
  87.      if (cudaStat != cudaSuccess) {
  88.       printf ("device memory allocation failed\n");
  89.       return EXIT_FAILURE;
  90.      }
  91.      cudaStat = cudaMalloc ((void**)&dev_D, k*l*sizeof(float));
  92.      if (cudaStat != cudaSuccess) {
  93.       printf ("device memory allocation failed\n");
  94.       return EXIT_FAILURE;
  95.      }
  96.      cudaStat = cudaMalloc ((void**)&dev_E, m*l*sizeof(float));
  97.      if (cudaStat != cudaSuccess) {
  98.       printf ("device memory allocation failed\n");
  99.       return EXIT_FAILURE;
  100.      }
  101.      
  102.      //copying matriсes
  103.      cudaMemcpy(dev_B, host_B, n*k*sizeof(float), cudaMemcpyHostToDevice);
  104.      cudaMemcpy(dev_A, host_A, m*n*sizeof(float), cudaMemcpyHostToDevice);
  105.      cudaMemcpy(dev_D, host_D, k*l*sizeof(float), cudaMemcpyHostToDevice);
  106.      float result[2];
  107.      
  108.      //GPU multiplication
  109.      //scalars
  110.      float scal[2];
  111.      scal[0]=1;
  112.      scal[1]=0;
  113.      
  114.      stat = cublasSgemm(handle,
  115.             CUBLAS_OP_N, CUBLAS_OP_N,
  116.             k, n, m,
  117.             scal,      // alpha (1)
  118.             dev_B, k,  // ?
  119.             dev_A, n,  // ?
  120.             (scal+1),  // beta (0)
  121.             dev_C, k); // ?
  122.      
  123.      if (stat != CUBLAS_STATUS_SUCCESS) {
  124.       printf ("Multiplication failed. (1)\n");
  125.       cudaFree (dev_A); cudaFree (dev_B);cudaFree (dev_C);
  126.       delete[] host_A; delete[] host_B; delete[] host_C;
  127.       cublasDestroy(handle);        
  128.       return EXIT_FAILURE;
  129.      }    
  130.      //CPU multiplication
  131.      mMul(m,n,k,host_A,host_B,host_C);
  132.  
  133.      
  134.      cudaThreadSynchronize();
  135.      
  136.      //copying result from device to host
  137.      host_result=new float[m*k];
  138.      stat = cublasGetVector(m*k, sizeof(float), dev_C, 1, host_result, 1);
  139.      if (stat != CUBLAS_STATUS_SUCCESS) {
  140.       printf ("data upload failed\n");
  141.       cudaFree (dev_A); cudaFree (dev_B); cudaFree (dev_C);
  142.       delete[] host_A; delete[] host_B; delete[] host_C;
  143.       delete[] host_result;
  144.       cublasDestroy(handle);        
  145.       return EXIT_FAILURE;
  146.      }    
  147.      
  148.      //comparing
  149.      
  150.      result[0]=0;
  151.      result[1]=0;
  152.      bool flag = false;
  153.      for (int i=0; i<m*k; i++)
  154.      {
  155.       result[0]+=host_result[i];
  156.       result[1]+=host_C[i];
  157.       if (abs(host_result[i] - host_C[i]) > 0.001)
  158.       {
  159.            flag=true;
  160.            break;
  161.       }
  162.      }
  163.      //if(result[0]!=result[1])
  164.      if (flag)
  165.       printf("Device: %f; Host: %f \n):\n", result[0], result[1]);
  166.      else
  167.       printf("(:\n");
  168.      
  169.      //The second multiplication
  170.      //CPU
  171.      stat = cublasSgemm(handle,
  172.             CUBLAS_OP_N, CUBLAS_OP_N,
  173.             l, k, m,
  174.                 scal,      // alpha (1)
  175.             dev_D, l,  // ?
  176.             dev_C, k,  // ?
  177.             (scal+1),  // beta (0)
  178.             dev_E, l); // ?
  179.      
  180.      if (stat != CUBLAS_STATUS_SUCCESS) {
  181.       printf ("Multiplication failed (2)\n");
  182.       cudaFree (dev_A);
  183.       cudaFree (dev_B);
  184.       cudaFree (dev_C);
  185.       delete[] host_A;
  186.       delete[] host_B;
  187.       delete[] host_C;
  188.       cublasDestroy(handle);        
  189.       return EXIT_FAILURE;
  190.      }    
  191.      
  192.      //CPU multiplication
  193.       mMul(m,k,l,host_C,host_D,host_E);
  194.      
  195.       //getting result
  196.       float *host_result1=new float[m*l];
  197.       cudaThreadSynchronize();
  198.  
  199.       //copying result to host
  200.       stat = cublasGetVector(m*l, sizeof(float), dev_E, 1, host_result1, 1);
  201.       if (stat != CUBLAS_STATUS_SUCCESS) {
  202.        printf ("data upload failed\n");
  203.        cudaFree (dev_A); cudaFree (dev_B); cudaFree (dev_C);
  204.        delete[] host_A; delete[] host_B; delete[] host_C;
  205.        delete[] host_result;
  206.        cublasDestroy(handle);        
  207.        return EXIT_FAILURE;
  208.       }    
  209.            
  210.       //comparing
  211.       result[0]=0;
  212.       result[1]=0;
  213.       for (int i=0; i<l*m; i++)
  214.       {
  215.        result[1]+=host_E[i];
  216.        result[0]+=host_result1[i];
  217.        if (abs(host_E[i] - host_result1[i]) > 0.001)
  218.        {
  219.         flag=true;
  220.            
  221.        }
  222.       }
  223.       //if(result[0]!=result[1])
  224.       if (flag)
  225.        printf("Device: %f; Host: %f \n):\n", result[0], result[1]);
  226.       else
  227.        printf("(:\n");
  228.      
  229.      
  230.       delete[] host_A;
  231.       delete[] host_B;
  232.       delete[] host_C;
  233.      
  234.       cudaFree (dev_A);
  235.       cudaFree (dev_B);
  236.       cudaFree (dev_C);
  237.      
  238.       return (EXIT_SUCCESS);
  239. }
Advertisement
RAW Paste Data Copied
Advertisement