title: Matrix inverse usng linear system solver through cublas, cublasCreate exception or something else
body: <p>I'm trying to inverse a matrix using linear equation solver through cublas CUDA library.</p>
  
  <p>The original equation is in form of:</p>
  
  <pre><code>Ax  = B  = I 
  
  I - identity matrix 
  A - The matrix I'm trying to inverse 
  x - (hopefully) the inverse(A) matrix
  </code></pre>
  
  <p>I'd like to perform LU decomposition, receive the following: </p>
  
  <pre><code>LUx = B 
  
  L - is a lower triangle matrix 
  U - is a upper triangle matrix 
  </code></pre>
  
  <p>Here is a good example who can show what I'm trying to do:</p>
  
  <p><a href="http://www.youtube.com/watch?v=dza5JTvMpzk" rel="nofollow noreferrer">http://www.youtube.com/watch?v=dza5JTvMpzk</a></p>
  
  <p>For the sake of discussion, <strong>let's assume I already have LU decomposition</strong> of A. (A = LU). I'm trying to find the inverse  in a two steps series:</p>
  
  <pre><code>Ax = B = I 
  
  LUx = B = I 
  
  1st step:  Declare y: 
  
  **Ly  = I** 
  
  solve 1st linear equation 
  
  2nd step, Solve the following linear equation
  
  **Ux = y** 
  
  find X = inverse(A) - *Bingo!@!* 
  </code></pre>
  
  <p>For now I have no idea what am I doing wrong here. There might be two assumptions, either I'm not using cublas properly or cublas throws an exception for no reason.</p>
  
  <p>See my code attached, there is Matlab code at the end:</p>
  
  <pre><code>     #include "cuda_runtime.h" 
       #include "device_launch_parameters.h"
  
       #include &lt;stdio.h&gt;
  
  
       //#include "cublas.h"
       #include "cublas_v2.h"
  
  
   int main()
   {
  
  cudaError_t cudaStat;
  cublasStatus_t stat;
  cublasHandle_t handle;
  
  //cublasInit();
  
  stat = cublasCreate(&amp;handle);
  if (stat != CUBLAS_STATUS_SUCCESS)
  {
      printf ("CUBLAS initialization failed\n");
      return -1;
  }
  
  unsigned int size = 3; 
  // Allowcate device matrixes 
  float* p_h_LowerTriangle; 
  float* p_d_LowerTriangle; 
  
  p_h_LowerTriangle = new float[size*size];
  p_h_LowerTriangle[0] = 1.f;  
  p_h_LowerTriangle[1] = 0.f;
  p_h_LowerTriangle[2] = 0.f;
  p_h_LowerTriangle[3] = 2.56f;
  p_h_LowerTriangle[4] = 1.f; 
  p_h_LowerTriangle[5] = 0.f; 
  p_h_LowerTriangle[6] = 5.76f; 
  p_h_LowerTriangle[7] = 3.5f;
  p_h_LowerTriangle[8] = 1.f;
  
  float* p_h_UpperTriangle; 
  float* p_d_UpperTriangle; 
  
  p_h_UpperTriangle = new float[size*size];
  p_h_UpperTriangle[0] = 25.f;  
  p_h_UpperTriangle[1] = 5.f;
  p_h_UpperTriangle[2] = 1.f;
  p_h_UpperTriangle[3] = 0.f;
  p_h_UpperTriangle[4] = -4.8f; 
  p_h_UpperTriangle[5] = -1.56f; 
  p_h_UpperTriangle[6] = 0.f; 
  p_h_UpperTriangle[7] = 0.f;
  p_h_UpperTriangle[8] = 0.7f;
  
  float* p_h_IdentityMatrix; 
  float* p_d_IdentityMatrix; 
  
  p_h_IdentityMatrix = new float[size*size];
  p_h_IdentityMatrix[0] = 1.f;  
  p_h_IdentityMatrix[1] = 0.f;
  p_h_IdentityMatrix[2] = 0.f;
  p_h_IdentityMatrix[3] = 0.f;
  p_h_IdentityMatrix[4] = 1.f; 
  p_h_IdentityMatrix[5] = 0.f; 
  p_h_IdentityMatrix[6] = 0.f; 
  p_h_IdentityMatrix[7] = 0.f;
  p_h_IdentityMatrix[8] = 1.f;
  
  cudaMalloc ((void**)&amp;p_d_LowerTriangle, size*size*sizeof(float));
  cudaMalloc ((void**)&amp;p_d_UpperTriangle, size*size*sizeof(float));
  cudaMalloc ((void**)&amp;p_d_IdentityMatrix, size*size*sizeof(float));
  
  
  float* p_d_tempMatrix; 
  cudaMalloc ((void**)&amp;p_d_tempMatrix, size*size*sizeof(float));
  
  
  stat = cublasSetMatrix(size,size,sizeof(float),p_h_LowerTriangle,size,p_d_LowerTriangle,size);
  stat = cublasSetMatrix(size,size,sizeof(float),p_h_UpperTriangle,size,p_d_UpperTriangle,size);
  stat = cublasSetMatrix(size,size,sizeof(float),p_h_IdentityMatrix,size,p_d_IdentityMatrix,size);
  
  cudaDeviceSynchronize();
  
  // 1st phase - solve Ly = b 
  const float alpha  = 1.f;
  
  // Function solves the triangulatr linear system with multiple right hand sides, function overrides b as a result 
  stat =  cublasStrsm(handle,
      cublasSideMode_t::CUBLAS_SIDE_LEFT, 
      cublasFillMode_t::CUBLAS_FILL_MODE_LOWER,
      cublasOperation_t::CUBLAS_OP_N,
      CUBLAS_DIAG_UNIT,
      size,
      size,
      &amp;alpha,
      p_d_LowerTriangle,
      size,
      p_d_IdentityMatrix,
      size);
  
  ////////////////////////////////////
  // TODO: printf 1st phase the results 
  cudaDeviceSynchronize();
  
  cudaMemcpy(p_h_IdentityMatrix,p_d_IdentityMatrix,size*size*sizeof(float),cudaMemcpyDeviceToHost);
  cudaMemcpy(p_h_UpperTriangle,p_d_UpperTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost);
  cudaMemcpy(p_h_LowerTriangle,p_d_LowerTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost);
  
  stat =cublasGetMatrix(size,size,sizeof(float),p_d_IdentityMatrix,size,p_h_IdentityMatrix,size);
  stat =cublasGetMatrix(size,size,sizeof(float),p_d_UpperTriangle,size,p_h_UpperTriangle,size);
  stat =cublasGetMatrix(size,size,sizeof(float),p_d_LowerTriangle,size,p_h_LowerTriangle,size);
  
  ////////////////////////////////////
  
  // 2nd phase - solve Ux = b
  stat =  cublasStrsm(handle,
      cublasSideMode_t::CUBLAS_SIDE_LEFT, 
      cublasFillMode_t::CUBLAS_FILL_MODE_UPPER,
      cublasOperation_t::CUBLAS_OP_N,
      CUBLAS_DIAG_NON_UNIT,
      size,
      size,
      &amp;alpha,
      p_d_UpperTriangle,
      size,
      p_d_IdentityMatrix,
      size);
  
  // TODO print the results 
  cudaDeviceSynchronize();
  
  ////////////////////////////////////
  cudaMemcpy(p_h_IdentityMatrix,p_d_IdentityMatrix,size*size*sizeof(float),cudaMemcpyDeviceToHost);
  cudaMemcpy(p_h_UpperTriangle,p_d_UpperTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost);
  cudaMemcpy(p_h_LowerTriangle,p_d_LowerTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost);
  ////////////////////////////////////
  
  cublasDestroy(handle);
  //cublasShutdown();
  
  cudaFree(p_d_LowerTriangle);
  cudaFree(p_d_UpperTriangle);
  cudaFree(p_d_IdentityMatrix);
  
  free(p_h_LowerTriangle);
  free(p_h_UpperTriangle);
  free(p_h_IdentityMatrix);
  
  
  return 0;
  }
  </code></pre>
  
  <p>Matlab code - works perfect:</p>
  
  <pre><code>  clear all
  
     UpperMatrix  = [ 25 5 1  ;  0 -4.8  -1.56  ;  0 0 0.7 ]
  
     LowerMatrix  = [1 0 0 ; 2.56 1 0 ; 5.76 3.5 1  ]
  
     IdentityMatrix = eye(3)
  
      % 1st phase solution 
      otps1.LT = true;
      y = linsolve(LowerMatrix,IdentityMatrix,otps1)
  
      %2nd phase solution 
      otps2.UT = true;
      y = linsolve(UpperMatrix,y,otps2)
  </code></pre>
  
  <h2>MATLAB output</h2>
  
  <p>UpperMatrix =</p>
  
  <pre><code>25.0000    5.0000    1.0000
        0   -4.8000   -1.5600
        0         0    0.7000
  </code></pre>
  
  <p>LowerMatrix =</p>
  
  <pre><code>1.0000         0         0
  2.5600    1.0000         0
  5.7600    3.5000    1.0000
  </code></pre>
  
  <p>IdentityMatrix =</p>
  
  <pre><code> 1     0     0
   0     1     0
   0     0     1
  </code></pre>
  
  <p>y =</p>
  
  <pre><code> 1.0000         0         0
  -2.5600    1.0000         0
   3.2000   -3.5000    1.0000
  </code></pre>
  
  <p>y =</p>
  
  <pre><code> 0.0476   -0.0833    0.0357
  -0.9524    1.4167   -0.4643
   4.5714   -5.0000    1.4286
  </code></pre>
  
  <p>UpperMatrix =</p>
  
  <pre><code>25.0000    5.0000    1.0000
        0   -4.8000   -1.5600
        0         0    0.7000
  </code></pre>
  
  <p>LowerMatrix =</p>
  
  <pre><code>1.0000         0         0
  2.5600    1.0000         0
  5.7600    3.5000    1.0000
  </code></pre>
  
  <p>IdentityMatrix =</p>
  
  <pre><code> 1     0     0
   0     1     0
   0     0     1
  </code></pre>
  
  <p>y =</p>
  
  <pre><code> 1.0000         0         0
  -2.5600    1.0000         0
   3.2000   -3.5000    1.0000
  
  &gt;&gt; 
  &gt;&gt; 
  &gt;&gt; 
  &gt;&gt; Inverse_LU_UT
  </code></pre>
  
  <p>UpperMatrix =</p>
  
  <pre><code>25.0000    5.0000    1.0000
        0   -4.8000   -1.5600
        0         0    0.7000
  </code></pre>
  
  <p>LowerMatrix =</p>
  
  <pre><code>1.0000         0         0
  2.5600    1.0000         0
  5.7600    3.5000    1.0000
  </code></pre>
  
  <p>IdentityMatrix =</p>
  
  <pre><code> 1     0     0
   0     1     0
   0     0     1
  </code></pre>
  
  <p>y =</p>
  
  <pre><code> 1.0000         0         0
  -2.5600    1.0000         0
   3.2000   -3.5000    1.0000
  </code></pre>
  
  <p>y =</p>
  
  <pre><code> 0.0476   -0.0833    0.0357
  -0.9524    1.4167   -0.4643
   4.5714   -5.0000    1.4286
  </code></pre>
  
  <p>I have no idea what am I doing wrong here, I suspect the cublasCreate  operation. Every time I hit the  command: </p>
  
  <pre><code>  cublasStatus_t stat;
    cublasHandle_t handle;
    stat = cublasCreate(&amp;handle); 
  </code></pre>
  
  <p>stat and handle variables are both valid. </p>
  
  <p>But VS10 output several error messages specified the following: </p>
  
  <blockquote>
    <p>First chance exception at 0x... microsoft C++ exception cudaError_enum at memory location 0x...</p>
  </blockquote>
  
  <p>Some may say it's an internal cublas error message, handled by the library itself. </p>