title: CUBLAS works unpredictably
body: <p>Wrote my first program using CUDA+CUBLAS. It just uses a 'cublasDgemm' function and computes a product of 2 N*N matrices.</p>
  
  <p>However, all the time I was launching my program, it keeped producing the same wrong answer (e.g. when multiplying 1*1 matrix containing 5 as a single element by 1*1 matrix containing element 6, it always said the result is 36, not 30).
  I checked the program several times with no success. But, when I came back to it the nexy day (i.e. after reboot), it worked just fine. I don't remember whether I recompiled it or not, but the truth is that it is the same VS project, same code, same computer with its GPU.</p>
  
  <p>So, can anyone explain me why could that have happened? And do I have to expect same strange behaviour further?</p>
  
  <p>Here is the code I was launching:</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;string&gt;
  #include &lt;iomanip&gt;
  #include &lt;cuda_runtime.h&gt;
  #include &lt;cublas_v2.h&gt;
  
  const int N = 5;
  #define IDX2F(i,j) ((i) * N + j)
  
  void fail(const cudaError_t&amp; cudaStatus, const std::string&amp; errorMessage) {
      if (cudaStatus != cudaSuccess) {
          std::cerr &lt;&lt; errorMessage &lt;&lt; std::endl;
          exit(EXIT_FAILURE);
      }
  }
  
  void fail(const cublasStatus_t&amp; status, const std::string&amp; errorMessage) {
      if (status != CUBLAS_STATUS_SUCCESS) {
          std::cerr &lt;&lt; errorMessage &lt;&lt; std::endl;
          exit(EXIT_FAILURE);
      }
  }
  
  void printMatrix(const double *C) {
      for (int i=0; i&lt;N; i++) {
          for (int j=0; j&lt;N; j++) {
              std::cout &lt;&lt; std::fixed &lt;&lt; std::setprecision(2) &lt;&lt; C[IDX2F(i,j)] &lt;&lt; ' ';
          }
          std::cout &lt;&lt; std::endl;
      }
      std::cout &lt;&lt; std::endl;
  }
  
  int main(int argc, char **argv) {
      cudaError_t cudaStatus;
      cublasStatus_t status;
      cublasHandle_t handle;
  
      double *A = new double[N*N];
      double *devPtrA;
  
      double *B = new double[N*N];
      double *devPtrB;
  
      double *C = new double[N*N];
      double *devPtrC;
  
      for (int i=0; i&lt;N; i++)
          for (int j=0; j&lt;N; j++)
              A[IDX2F(i,j)] = i + j;
  
      for (int i=0; i&lt;N; i++)
          for (int j=0; j&lt;N; j++)
              B[IDX2F(i,j)] = i + j * 0.5;
  
      // do not have to set anything into matrix C, because beta = 0
  
      // allocate mamory on GPU
      cudaStatus = cudaMalloc((void**)&amp;devPtrC, N*N*sizeof(*C));
      fail(cudaStatus, "device memory allocation failed");
  
      cudaStatus = cudaMalloc((void**)&amp;devPtrA, N*N*sizeof(*A));
      fail(cudaStatus, "device memory allocation failed");
  
      cudaStatus = cudaMalloc((void**)&amp;devPtrB, N*N*sizeof(*B));
      fail(cudaStatus, "device memory allocation failed");
  
      // create GPU handle
      status = cublasCreate(&amp;handle);
      fail(status, "CUBLAS initialization failed");
  
      // copying matrices from host to GPU
      status = cublasSetMatrix(N, N, sizeof (*B), B, N, devPtrB, N);
      fail(status, "failed to load data from host to GPU");
  
      status = cublasSetMatrix(N, N, sizeof (*A), A, N, devPtrA, N);
      fail(status, "failed to load data from host to GPU");
  
      const double ONE = 1;
      const double ZERO = 0;
  
      printMatrix(A);
      printMatrix(B);
  
      status = cublasDgemm(   handle,
                              CUBLAS_OP_N, CUBLAS_OP_N,
                              N, N, N,
                              &amp;ONE,
                              devPtrA, N,
                              devPtrB, N,
                              &amp;ZERO,
                              devPtrC, N);
  
      fail(status, "error cublasDgemm");
  
      status = cublasGetMatrix(N, N, sizeof (*C), devPtrC, N, C, N);
      fail(status, "could not load result back from GPU to host");
  
      printMatrix(C);
  
      status = cublasDestroy(handle);
      fail(status, "could not destroy CUBLAS handle");
  
      cudaStatus = cudaFree(devPtrC);
      fail(cudaStatus, "device memory freeing failed");
  
      cudaStatus = cudaFree(devPtrB);
      fail(cudaStatus, "device memory freeing failed");
  
      cudaStatus = cudaFree(devPtrA);
      fail(cudaStatus, "device memory freeing failed");
  
      delete[] C;
      delete[] B;
      delete[] A;
  
      return EXIT_SUCCESS;
  }
  </code></pre>