title: CUDA cuSolver gesvdj with large matrix
body: <p>I am running a slightly modified version of the code found <a href="https://docs.nvidia.com/cuda/cusolver/index.html" rel="nofollow noreferrer">here</a> under section "G.2. SVD with singular vectors (via Jacobi method)" on a NVIDIA P6000. The slight modifications are to dynamically allocate memory in the heap for the A, U, and V vectors and to fill the A vector of a specified size with values that depend on the index into A. I also converted everything from doubles to floats. The last modification is a loop over the gesvdj call itself and convergence checking up to some # of iterations (10 in my case). </p>
  
  <p>With those slight modifications, I was able to overcome my first hurdle of performing an SVD on symmetric arrays bigger than ~1000x1000 in size. I eventually need to run an SVD on a array of size 1048576x20. </p>
  
  <p>Currently, the algorithm runs for arrays of size 10000x20 but fails when I go to 50000x20. </p>
  
  <p>The issue seems to be stemming from the gesvdj call itself. A synchronize call after calling gesvdj is failing and returning a generic access error. </p>
  
  <p>If I run the program using cuda-memcheck I get a series of these errors for different threads in the same block:</p>
  
  <pre><code>Invalid __global__ write of size 4
  =========     at 0x00000108 in void eye_kernel&lt;float, int=5, int=3&gt;(int, int, float*, int)
  =========     by thread (0,7,0) in block (16,1342,0)
  =========     Address 0x7f4ed40414c0 is out of bounds
  =========     Saved host backtrace up to driver entry point at kernel launch time
  =========     Host Frame:/lib64/libcuda.so.1 (cuLaunchKernel + 0x2c5) [0x269e85]
  =========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x100c822]
  =========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x100ca17]
  =========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x1040dd5]
  =========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 [0x235c5d]
  =========     Host Frame:/usr/local/cuda-10.0/lib64/libcusolver.so.10.0 (cusolverDnSgesvdj + 0x508) [0x21f9a8]
  =========     Host Frame:./gesvdj_example [0x4518]
  =========     Host Frame:/lib64/libc.so.6 (__libc_start_main + 0xf5) [0x223d5]
  =========     Host Frame:./gesvdj_example [0x3ab9]
  </code></pre>
  
  <p>I'm wondering if I have hit some sort of cusolver internal limitation? Anyone have any ideas? I can provide my exact code if necessary, but its so similar to the example that I figured I would just point people there. </p>
  
  <p>Thanks!</p>
  
  <p>Edit to add the offending code from the example I linked to, the algorithm is failing at the assert(CUSOLVER_STATUS_SUCCESS == status); line. I'm really new to coding in C and CUDA so sorry if there is some obvious debugging information I am leaving out. </p>
  
  <pre><code>/* step 5: compute SVD */
      status = cusolverDnDgesvdj(
          cusolverH,
          jobz,  /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
                 /* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
          econ,  /* econ = 1 for economy size */
          m,     /* nubmer of rows of A, 0 &lt;= m */
          n,     /* number of columns of A, 0 &lt;= n  */
          d_A,   /* m-by-n */
          lda,   /* leading dimension of A */
          d_S,   /* min(m,n)  */
                 /* the singular values in descending order */
          d_U,   /* m-by-m if econ = 0 */
                 /* m-by-min(m,n) if econ = 1 */
          lda,   /* leading dimension of U, ldu &gt;= max(1,m) */
          d_V,   /* n-by-n if econ = 0  */
                 /* n-by-min(m,n) if econ = 1  */
          lda,   /* leading dimension of V, ldv &gt;= max(1,n) */
          d_work,
          lwork,
          d_info,
          gesvdj_params);
      cudaStat1 = cudaDeviceSynchronize();
      assert(CUSOLVER_STATUS_SUCCESS == status);
      assert(cudaSuccess == cudaStat1);
  </code></pre>
  
  <p><strong>Edit 2</strong> Adding my code...</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;string.h&gt;
  #include &lt;assert.h&gt;
  #include &lt;cusolverDn.h&gt;
  #include&lt;iostream&gt;
  #include&lt;iomanip&gt;
  #include&lt;assert.h&gt;
  #include &lt;cuda_runtime_api.h&gt;
  #include "cuda_runtime.h"
  
  int main(int argc, char*argv[])
  {
  
      //Setting which device to run on
      const int device = 0;
      cudaSetDevice(device);
  
      cusolverDnHandle_t cusolverH = NULL;
      cudaStream_t stream = NULL;
      gesvdjInfo_t gesvdj_params = NULL;
  
      cusolverStatus_t status = CUSOLVER_STATUS_SUCCESS;
      cudaError_t cudaStat1 = cudaSuccess;
      cudaError_t cudaStat2 = cudaSuccess;
      cudaError_t cudaStat3 = cudaSuccess;
      cudaError_t cudaStat4 = cudaSuccess;
      cudaError_t cudaStat5 = cudaSuccess;
      const long int m = 50000;
      const long int n = 20;
      const int lda = m;
  
      // --- Setting the host, Nrows x Ncols matrix
      float *A = (float *)malloc(m * n * sizeof(float));
        for(long int j = 0; j &lt; m; j++)
            for(long int i = 0; i &lt; n; i++)
                A[j + i*m] = sqrt((float)(i + j));
  
      float *U = (float *)malloc(m * m * sizeof(float)); /* m-by-m unitary matrix, left singular vectors  */
      float *V = (float *)malloc(m * n * sizeof(float)); /* n-by-n unitary matrix, right singular vectors */
      float S[n];     /* numerical singular value */ 
      float *d_A = NULL;  /* device copy of A */
      float *d_S = NULL;  /* singular values */
      float *d_U = NULL;  /* left singular vectors */
      float *d_V = NULL;  /* right singular vectors */
      int *d_info = NULL;  /* error info */
      int lwork = 0;       /* size of workspace */
      float *d_work = NULL; /* devie workspace for gesvdj */
      int info = 0;        /* host copy of error info */
  
  /* configuration of gesvdj  */
      const double tol = 1.e-7;
      const int max_sweeps = 100;
      const cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvectors.
      const int econ = 0 ; /* econ = 1 for economy size */
  
  /* numerical results of gesvdj  */
      double residual = 0;
      int executed_sweeps = 0;
  
      printf("tol = %E, default value is machine zero \n", tol);
      printf("max. sweeps = %d, default value is 100\n", max_sweeps);
      printf("econ = %d \n", econ);
      printf("=====\n");
  
  /* step 1: create cusolver handle, bind a stream */
      status = cusolverDnCreate(&amp;cusolverH);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
      cudaStat1 = cudaStreamCreateWithFlags(&amp;stream, cudaStreamNonBlocking);
      assert(cudaSuccess == cudaStat1);
  
      status = cusolverDnSetStream(cusolverH, stream);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
  /* step 2: configuration of gesvdj */
      status = cusolverDnCreateGesvdjInfo(&amp;gesvdj_params);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
  /* default value of tolerance is machine zero */
      status = cusolverDnXgesvdjSetTolerance(
          gesvdj_params,
          tol);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
  /* default value of max. sweeps is 100 */
      status = cusolverDnXgesvdjSetMaxSweeps(
          gesvdj_params,
          max_sweeps);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
  /* step 3: copy A to device */
      cudaStat1 = cudaMalloc ((void**)&amp;d_A   , sizeof(float)*lda*n);
      cudaStat2 = cudaMalloc ((void**)&amp;d_S   , sizeof(float)*n);
      cudaStat3 = cudaMalloc ((void**)&amp;d_U   , sizeof(float)*lda*m);
      cudaStat4 = cudaMalloc ((void**)&amp;d_V   , sizeof(float)*lda*n);
      cudaStat5 = cudaMalloc ((void**)&amp;d_info, sizeof(int));
      assert(cudaSuccess == cudaStat1);
      assert(cudaSuccess == cudaStat2);
      assert(cudaSuccess == cudaStat3);
      assert(cudaSuccess == cudaStat4);
      assert(cudaSuccess == cudaStat5);
  
      cudaStat1 = cudaMemcpy(d_A, A, sizeof(float)*lda*n, cudaMemcpyHostToDevice);
      assert(cudaSuccess == cudaStat1);
      /* step 4: query workspace of SVD */
      status = cusolverDnSgesvdj_bufferSize(
          cusolverH,
          jobz, /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
                /* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
          econ, /* econ = 1 for economy size */
          m,    /* nubmer of rows of A, 0 &lt;= m */
          n,    /* number of columns of A, 0 &lt;= n  */
          d_A,  /* m-by-n */
          lda,  /* leading dimension of A */
          d_S,  /* min(m,n) */
                /* the singular values in descending order */
          d_U,  /* m-by-m if econ = 0 */
                /* m-by-min(m,n) if econ = 1 */
          lda,  /* leading dimension of U, ldu &gt;= max(1,m) */
          d_V,  /* n-by-n if econ = 0  */
                /* n-by-min(m,n) if econ = 1  */
          lda,  /* leading dimension of V, ldv &gt;= max(1,n) */
          &amp;lwork,
          gesvdj_params);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
      cudaStat1 = cudaMalloc((void**)&amp;d_work , sizeof(float)*lwork);
      assert(cudaSuccess == cudaStat1);
  
  /* step 5: compute SVD */
      //Iterating over SVD calculation, not part of example
      int iters;
      for (iters = 10; iters &gt; 0; iters--){
  
      status = cusolverDnSgesvdj(
          cusolverH,
          jobz,  /* CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only */
                 /* CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors */
          econ,  /* econ = 1 for economy size */
          m,     /* nubmer of rows of A, 0 &lt;= m */
          n,     /* number of columns of A, 0 &lt;= n  */
          d_A,   /* m-by-n */
          lda,   /* leading dimension of A */
          d_S,   /* min(m,n)  */
                 /* the singular values in descending order */
          d_U,   /* m-by-m if econ = 0 */
                 /* m-by-min(m,n) if econ = 1 */
          lda,   /* leading dimension of U, ldu &gt;= max(1,m) */
          d_V,   /* n-by-n if econ = 0  */
                 /* n-by-min(m,n) if econ = 1  */
          lda,   /* leading dimension of V, ldv &gt;= max(1,n) */
          d_work,
          lwork,
          d_info,
          gesvdj_params);
  
      cudaStat1 = cudaDeviceSynchronize();
      assert(CUSOLVER_STATUS_SUCCESS == status);
      assert(cudaSuccess == cudaStat1);
  
      cudaStat1 = cudaMemcpy(U, d_U, sizeof(float)*lda*m, cudaMemcpyDeviceToHost);
      cudaStat2 = cudaMemcpy(V, d_V, sizeof(float)*lda*n, cudaMemcpyDeviceToHost);
      cudaStat3 = cudaMemcpy(S, d_S, sizeof(float)*n    , cudaMemcpyDeviceToHost);
      cudaStat4 = cudaMemcpy(&amp;info, d_info, sizeof(int), cudaMemcpyDeviceToHost);
      cudaStat5 = cudaDeviceSynchronize();
      assert(cudaSuccess == cudaStat1);
      assert(cudaSuccess == cudaStat2);
      assert(cudaSuccess == cudaStat3);
      assert(cudaSuccess == cudaStat4);
      assert(cudaSuccess == cudaStat5);
  
      if ( 0 == info ){
          printf("gesvdj converges \n");
      }else if ( 0 &gt; info ){
          printf("%d-th parameter is wrong \n", -info);
          exit(1);
      }else{
          printf("WARNING: info = %d : gesvdj does not converge \n", info );
      }
      }
  /* step 6: measure error of singular value */
      status = cusolverDnXgesvdjGetSweeps(
          cusolverH,
          gesvdj_params,
          &amp;executed_sweeps);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
      status = cusolverDnXgesvdjGetResidual(
          cusolverH,
          gesvdj_params,
          &amp;residual);
      assert(CUSOLVER_STATUS_SUCCESS == status);
  
      printf("residual |A - U*S*V**H|_F = %E \n", residual );
      printf("number of executed sweeps = %d \n\n", executed_sweeps );
  
  /*  free resources  */
      if (A      ) free(A);
      if (V      ) free(V);
      if (U      ) free(U);
      if (d_A    ) cudaFree(d_A);
      if (d_S    ) cudaFree(d_S);
      if (d_U    ) cudaFree(d_U);
      if (d_V    ) cudaFree(d_V);
      if (d_info ) cudaFree(d_info);
      if (d_work ) cudaFree(d_work);
  
      if (cusolverH    ) cusolverDnDestroy(cusolverH);
      if (stream       ) cudaStreamDestroy(stream);
      if (gesvdj_params) cusolverDnDestroyGesvdjInfo(gesvdj_params);
  
      cudaDeviceReset();
      return 0;
  }
  </code></pre>