title: Improving the solution of sparse linear systems
body: <p>I have written a code in C++ on a Linux system which solves the linear system <code>A x = b</code>, where <code>A</code> is a sparse symmetric matrix using the following two approaches:</p>
  
  <ol>
  <li>Using <code>UMFPACK</code> to sequentially factorize and do backward forward substitution.</li>
  <li>Using <code>UMFPACK</code> to sequentially factorize followed by a backward forward substitution using the <code>cuSPARSE</code> library.</li>
  </ol>
  
  <p>The system configuration I have is: CUDA version 5.0, <code>UMFPACK</code> version 5.6.2, Linux kernel version  Debian 3.2.46-1, graphics card used: GeForce GTX Titan.  </p>
  
  <p>Theoretically, the second approach should perform better than the first approach with minimal or no errors. However, I am observing the following problems:</p>
  
  <ol>
  <li>The backward/forward substitution using <code>UMFPACK</code> function <strong>umfpack_di_solve</strong> is almost <code>2x</code> faster than the CUDA variant.</li>
  <li>For some matrices, the error between the results obtained using <code>UMFPACK</code> and CUDA are quite large with the maximum error being <code>3.2537</code>, whereas, for other matrices, it is of the order of <code>1e-16</code>.</li>
  </ol>
  
  <p>Attached is my tar file with the following components:</p>
  
  <ul>
  <li>A folder factorize_copy with the main file <strong>fc.cu</strong> that I am using to solve the linear system. It reads the sparse matrices from the <strong>grid_*_CSC.m</strong> files which are also present in the same directory. For convenience, the results with the three sparse matrices provided are also given in the text files.  </li>
  <li>A folder with all the dependencies for compiling and running <code>UMFPACK</code> (which we are also using for computations). </li>
  </ul>
  
  <p>The link to the tar file is 
  <a href="https://www.dropbox.com/s/9qfs5awclshyk3b/code.tar.gz" rel="nofollow">https://www.dropbox.com/s/9qfs5awclshyk3b/code.tar.gz</a></p>
  
  <p>If you wish to run the code, I have provided a <strong>MAKEFILE</strong> as I use in my system in the factorize_copy directory. You may need to recompile the <code>UMFPACK</code> Libraries. </p>
  
  <p>A sample output from our program for a <code>586 x 586</code> Sparse Matrix is also shown below (Note that errors for this case are very high as compared to other sparse matrices that we checked).</p>
  
  <pre>
  ***** Reading the Grids
  
      Reading Grids Successful
  
  ***** Solving a sparse matrix of size: 586x586
  
  ***** Solving the grid on umfpack
  
  ***** Factorizing The Grid
  
  -------------- CPU TIME for umfpack factorization is: 0.00109107
  
  -------------- Wall-Clock TIME for umfpack factorization is: 0
  
      Factorizing the Grid successful
  
      Solving the grid on umfpack successful
  
  -------------- <b>CPU TIME for umfpack solve is: 6.281e-05</b>
  
  ***** Allocating GPU memory and Copying data
  
  ---------------- CPU Time for Allocating GPU memory and Copying Data: 1.6
  
  ***** Performing b = P*b to account for the row ordering in A 
  
      Matrix-Vector (Pb) multiplication successful
  
  ***** Solving the system: LUx=b
  
      Analyzing Ly = b successful
  
      Solving Ly = b successful
  
      Analyzing Ux = y successful
  
      Solving Ux = y successful
  
  ***** Performing x = Q*x to account for the column ordering in A
  
      Matrix-Vector (Qx) multiplication successful
  
  ---------- <b>GPU Time for the solve is: 5.68029 ms</b>
  
  ##### <b>Maximum error between UMFPACK and CUDA: 3.2537
  
  ##### Average error between UMFPACK and CUDA: 0.699926</b>
  
  ***** Writing the results to the output files
  
      Result Written to the file 'vout_586.m' and the file 'vout_umfpack_586.m'
  
  (Operation Successful!)
  
  </pre>
  
  <p>I would really appreciate if anyone could point out what the possible error could be in this case. If there is a better method of solving sparse linear systems using CUDA that I am missing out, please let me know. </p>
  
  <p><strong>EDIT:</strong> I figured out why it was giving error in some cases and in some cases not. I had a mistake in the number of threads per block when calling a kernel function in the code. However, I still have the issue of gaining a speed-up. </p>