title: Best way of solving sparse linear systems in C++ - GPU Possible?
body: <p>I am currently working on a project where we need to solve</p>
  
  <p><code>|Ax - b|^2</code>.</p>
  
  <p>In this case, <code>A</code> is a very sparse matrix and <code>A'A</code> has at most 5 nonzero elements in each row. </p>
  
  <p>We are working with images and the dimension of <code>A'A</code> is <code>NxN</code> where N is the number of pixels. In this case <code>N = 76800</code>. We plan to go to <code>RGB</code> and then the dimension will be <code>3Nx3N</code>. </p>
  
  <p>In matlab solving <code>(A'A)\(A'b)</code> takes about 0.15 s, using doubles. </p>
  
  <p>I have now done some experimenting with <code>Eigens</code> sparse solvers. I have tried:</p>
  
  <pre><code>SimplicialLLT
  SimplicialLDLT
  SparseQR
  ConjugateGradient
  </code></pre>
  
  <p>and some different orderings. The by far best so far is </p>
  
  <pre><code>SimplicialLDLT
  </code></pre>
  
  <p>which takes about <code>0.35 - 0.5</code> using <code>AMDOrdering</code>. </p>
  
  <p>When I for example use <code>ConjugateGradient</code> it takes roughly <code>6 s</code>, using <code>0</code> as initilization.</p>
  
  <p>The code for solving the problem is:    </p>
  
  <pre><code>    A_tot.makeCompressed();
       // Create solver
      Eigen::SimplicialLDLT&lt;Eigen::SparseMatrix&lt;float&gt;, Eigen::Lower, Eigen::AMDOrdering&lt;int&gt; &gt; solver;
      // Eigen::ConjugateGradient&lt;Eigen::SparseMatrix&lt;float&gt;, Eigen::Lower&gt; cg;
      solver.analyzePattern(A_tot);
      t1 = omp_get_wtime();
      solver.compute(A_tot);
      if (solver.info() != Eigen::Success)
       {
           std::cerr &lt;&lt; "Decomposition Failed" &lt;&lt; std::endl;
           getchar();
       }
      Eigen::VectorXf opt = solver.solve(b_tot);
      t2 = omp_get_wtime();
      std::cout &lt;&lt; "Time for normal equations: " &lt;&lt; t2 - t1 &lt;&lt; std::endl;
  </code></pre>
  
  <p>This is the first time I use sparse matrices in C++ and its solvers. For this project speed is crucial and below <code>0.1 s</code> is a minimum. </p>
  
  <p>I would like to get some feedback on which would be the best strategy here. For example one is supposed to be able to use <code>SuiteSparse</code> and <code>OpenMP</code> in <code>Eigen</code>. What are your experiences about these types of problems? Is there a way of extracting the structure for example? And should <code>conjugateGradient</code> really be that slow?</p>
  
  <p>Edit:</p>
  
  <p>Thanks for som valuable comments! Tonight I have been reading a bit about cuSparse on Nvidia. It seems to be able to do factorisation an even solve systems. In particular it seems one can reuse pattern and so forth. The question is how fast could this be and what is the possible overhead?</p>
  
  <p>Given that the amount of data in my matrix A is the same as in an image, the memory copying should not be such an issue. I did some years ago software for real-time 3D reconstruction and then you copy data for each frame and a slow version still runs in over 50 Hz. So if the factorization is much faster it is a possible speed-up? In particualar the rest of the project will be on the GPU, so if one can solve it there directly and keep the solution it is no drawback I guess.</p>
  
  <p>A lot has happened in the field of Cuda and I am not really up to date.</p>
  
  <p>Here are two links I found: <a href="https://developer.nvidia.com/cusolver" rel="nofollow noreferrer">Benchmark?</a>, <a href="https://se.mathworks.com/help/distcomp/examples/benchmarking-a-b-on-the-gpu.html" rel="nofollow noreferrer">MatlabGPU</a></p>