title: Solving sparse definite positive linear systems in CUDA
body: <p>We are experiencing problems while using <code>cuSOLVER</code>'s <code>cusolverSpScsrlsvchol</code> function, probably due to misunderstanding of the <code>cuSOLVER</code> library.</p>
  
  <p>Motivation: we are solving the Poisson equation <code>-divgrad x = b</code> on a rectangular grid. In <code>2</code> dimensions with a <code>5</code>-stencil <code>(1, 1, -4, 1, 1)</code>, the Laplacian on the grid provides a (quite sparse) matrix <code>A</code>. Moreover, the charge distribution on the grid gives a (dense) vector <code>b</code>. <code>A</code> is positive definite and symmetric. </p>
  
  <p>Now we solve <code>A * x = b</code> for <code>x</code> using nvidia's new <code>cuSOLVER</code> library that comes with CUDA 7.0 . It provides a function <code>cusolverSpScsrlsvchol</code> which should do the sparse Cholesky factorisation for floats.</p>
  
  <p>Note: we are able to correctly solve the system with the alternative sparse QR factorisation function <code>cusolverSpScsrlsvqr</code>. For a <code>4 x 4</code> grid with all <code>b</code> entries on the edge being <code>1</code> and the rest <code>0</code>, we get for <code>x</code>:</p>
  
  <pre><code>1 1 0.999999 1 1 1 0.999999 1 1 1 1 1 1 1 1 1 
  </code></pre>
  
  <p>Our problems:</p>
  
  <ol>
  <li><p><code>cusolverSpScsrlsvchol</code> returns wrong results for <code>x</code>:</p>
  
  <pre><code>1 3.33333 2.33333 1 3.33333 2.33333 1.33333 1 2.33333 1.33333 0.666667 1 1 1 1 1
  </code></pre></li>
  <li><p>(solved, see answer below) Converting the CSR matrix <code>A</code> to a dense matrix and showing the output gives weird numbers (<code>10^-44</code> and the like). The respective data from the CSR format are correct and validated with python numpy.</p></li>
  <li><p>(solved, see answer below) The alternative sparse <code>LU</code> and partial pivoting with <code>cusolverSpScsrlsvlu</code> cannot even be found:</p>
  
  <pre><code>$ nvcc -std=c++11 cusparse_test3.cu -o cusparse_test3 -lcusparse -lcusolver
  cusparse_test3.cu(208): error: identifier "cusolverSpScsrlsvlu" is undefined
  </code></pre></li>
  </ol>
  
  <p>What are we doing wrong? Thanks for your help!</p>
  
  <p>Our C++ CUDA code:</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;cuda_runtime.h&gt;
  #include &lt;cuda.h&gt;
  #include &lt;cusolverSp.h&gt;
  #include &lt;cusparse.h&gt;
  #include &lt;vector&gt;
  #include &lt;cassert&gt;
  
  
  // create poisson matrix with Dirichlet bc. of a rectangular grid with
  // dimension NxN
  void assemble_poisson_matrix_coo(std::vector&lt;float&gt;&amp; vals, std::vector&lt;int&gt;&amp; row, std::vector&lt;int&gt;&amp; col,
                       std::vector&lt;float&gt;&amp; rhs, int Nrows, int Ncols) {
  
          //nnz: 5 entries per row (node) for nodes in the interior
      // 1 entry per row (node) for nodes on the boundary, since we set them explicitly to 1.
      int nnz = 5*Nrows*Ncols - (2*(Ncols-1) + 2*(Nrows-1))*4;
      vals.resize(nnz);
      row.resize(nnz);
      col.resize(nnz);
      rhs.resize(Nrows*Ncols);
  
      int counter = 0;
      for(int i = 0; i &lt; Nrows; ++i) {
          for (int j = 0; j &lt; Ncols; ++j) {
              int idx = j + Ncols*i;
              if (i == 0 || j == 0 || j == Ncols-1 || i == Nrows-1) {
                  vals[counter] = 1.;
                  row[counter] = idx;
                  col[counter] = idx;
                  counter++;
                  rhs[idx] = 1.;
  //                if (i == 0) {
  //                    rhs[idx] = 3.;
  //                }
              } else { // -laplace stencil
                  // above
                  vals[counter] = -1.;
                  row[counter] = idx;
                  col[counter] = idx-Ncols;
                  counter++;
                  // left
                  vals[counter] = -1.;
                  row[counter] = idx;
                  col[counter] = idx-1;
                  counter++;
                  // center
                  vals[counter] = 4.;
                  row[counter] = idx;
                  col[counter] = idx;
                  counter++;
                  // right
                  vals[counter] = -1.;
                  row[counter] = idx;
                  col[counter] = idx+1;
                  counter++;
                  // below
                  vals[counter] = -1.;
                  row[counter] = idx;
                  col[counter] = idx+Ncols;
                  counter++;
  
                  rhs[idx] = 0;
              }
          }
      }
      assert(counter == nnz);
  }
  
  
  
  int main() {
      // --- create library handles:
      cusolverSpHandle_t cusolver_handle;
      cusolverStatus_t cusolver_status;
      cusolver_status = cusolverSpCreate(&amp;cusolver_handle);
      std::cout &lt;&lt; "status create cusolver handle: " &lt;&lt; cusolver_status &lt;&lt; std::endl;
  
      cusparseHandle_t cusparse_handle;
      cusparseStatus_t cusparse_status;
      cusparse_status = cusparseCreate(&amp;cusparse_handle);
      std::cout &lt;&lt; "status create cusparse handle: " &lt;&lt; cusparse_status &lt;&lt; std::endl;
  
      // --- prepare matrix:
      int Nrows = 4;
      int Ncols = 4;
      std::vector&lt;float&gt; csrVal;
      std::vector&lt;int&gt; cooRow;
      std::vector&lt;int&gt; csrColInd;
      std::vector&lt;float&gt; b;
  
      assemble_poisson_matrix_coo(csrVal, cooRow, csrColInd, b, Nrows, Ncols);
  
      int nnz = csrVal.size();
      int m = Nrows * Ncols;
      std::vector&lt;int&gt; csrRowPtr(m+1);
  
      // --- prepare solving and copy to GPU:
      std::vector&lt;float&gt; x(m);
      float tol = 1e-5;
      int reorder = 0;
      int singularity = 0;
  
      float *db, *dcsrVal, *dx;
      int *dcsrColInd, *dcsrRowPtr, *dcooRow;
      cudaMalloc((void**)&amp;db, m*sizeof(float));
      cudaMalloc((void**)&amp;dx, m*sizeof(float));
      cudaMalloc((void**)&amp;dcsrVal, nnz*sizeof(float));
      cudaMalloc((void**)&amp;dcsrColInd, nnz*sizeof(int));
      cudaMalloc((void**)&amp;dcsrRowPtr, (m+1)*sizeof(int));
      cudaMalloc((void**)&amp;dcooRow, nnz*sizeof(int));
  
      cudaMemcpy(db, b.data(), b.size()*sizeof(float), cudaMemcpyHostToDevice);
      cudaMemcpy(dcsrVal, csrVal.data(), csrVal.size()*sizeof(float), cudaMemcpyHostToDevice);
      cudaMemcpy(dcsrColInd, csrColInd.data(), csrColInd.size()*sizeof(int), cudaMemcpyHostToDevice);
      cudaMemcpy(dcooRow, cooRow.data(), cooRow.size()*sizeof(int), cudaMemcpyHostToDevice);
  
      cusparse_status = cusparseXcoo2csr(cusparse_handle, dcooRow, nnz, m,
                                         dcsrRowPtr, CUSPARSE_INDEX_BASE_ZERO);
      std::cout &lt;&lt; "status cusparse coo2csr conversion: " &lt;&lt; cusparse_status &lt;&lt; std::endl;
  
      cudaDeviceSynchronize(); // matrix format conversion has to be finished!
  
      // --- everything ready for computation:
  
      cusparseMatDescr_t descrA;
  
      cusparse_status = cusparseCreateMatDescr(&amp;descrA);
      std::cout &lt;&lt; "status cusparse createMatDescr: " &lt;&lt; cusparse_status &lt;&lt; std::endl;
  
      // optional: print dense matrix that has been allocated on GPU
  
      std::vector&lt;float&gt; A(m*m, 0);
      float *dA;
      cudaMalloc((void**)&amp;dA, A.size()*sizeof(float));
  
      cusparseScsr2dense(cusparse_handle, m, m, descrA, dcsrVal,
                         dcsrRowPtr, dcsrColInd, dA, m);
  
      cudaMemcpy(A.data(), dA, A.size()*sizeof(float), cudaMemcpyDeviceToHost);
      std::cout &lt;&lt; "A: \n";
      for (int i = 0; i &lt; m; ++i) {
          for (int j = 0; j &lt; m; ++j) {
              std::cout &lt;&lt; A[i*m + j] &lt;&lt; " ";
          }
          std::cout &lt;&lt; std::endl;
      }
  
      cudaFree(dA);
  
      std::cout &lt;&lt; "b: \n";
      cudaMemcpy(b.data(), db, (m)*sizeof(int), cudaMemcpyDeviceToHost);
      for (auto a : b) {
          std::cout &lt;&lt; a &lt;&lt; ",";
      }
      std::cout &lt;&lt; std::endl;
  
  
      // --- solving!!!!
  
  //    cusolver_status = cusolverSpScsrlsvchol(cusolver_handle, m, nnz, descrA, dcsrVal,
  //                       dcsrRowPtr, dcsrColInd, db, tol, reorder, dx,
  //                       &amp;singularity);
  
       cusolver_status = cusolverSpScsrlsvqr(cusolver_handle, m, nnz, descrA, dcsrVal,
                          dcsrRowPtr, dcsrColInd, db, tol, reorder, dx,
                          &amp;singularity);
  
      cudaDeviceSynchronize();
  
      std::cout &lt;&lt; "singularity (should be -1): " &lt;&lt; singularity &lt;&lt; std::endl;
  
      std::cout &lt;&lt; "status cusolver solving (!): " &lt;&lt; cusolver_status &lt;&lt; std::endl;
  
      cudaMemcpy(x.data(), dx, m*sizeof(float), cudaMemcpyDeviceToHost);
  
      // relocated these 2 lines from above to solve (2):
      cusparse_status = cusparseDestroy(cusparse_handle);
      std::cout &lt;&lt; "status destroy cusparse handle: " &lt;&lt; cusparse_status &lt;&lt; std::endl;
  
      cusolver_status = cusolverSpDestroy(cusolver_handle);
      std::cout &lt;&lt; "status destroy cusolver handle: " &lt;&lt; cusolver_status &lt;&lt; std::endl;
  
      for (auto a : x) {
          std::cout &lt;&lt; a &lt;&lt; " ";
      }
      std::cout &lt;&lt; std::endl;
  
  
  
      cudaFree(db);
      cudaFree(dx);
      cudaFree(dcsrVal);
      cudaFree(dcsrColInd);
      cudaFree(dcsrRowPtr);
      cudaFree(dcooRow);
  
      return 0;
  }
  </code></pre>