title: Segmentation fault (core dumped) with my mex cuda code while exiting MATLAB
body: <p>I have a mixed mex and cuda code to evaluate phi = 1/2*x'*A*x - b'*x, where x and b are both m by 1 vectors and A is an m by m matrix. The code can be compiled and executed and it also gives me the correct answer. </p>
  
  <p>However, when I exited MATLAB, I kept getting the error Segmentation fault (core dumped). What I did in the code is that, I generate A, b and x in MATLAB, use a mex function to pass them to cuda. Then I evaluate phi = 1/2*x'*A*x - b'*x on GPU (use the cuda linear algebra library cublas) and use mex to transfer phi back to MATLAB. </p>
  
  <p>Can anyone help me to see where the problem is? Thanks in advance.</p>
  
  <p>Btw, here is how I compile it:</p>
  
  <pre><code>nvcc -arch=sm_20 -c test.cu -Xcompiler -fPIC -I/site/local/matlab-r2012a/extern/include/
  mex -L/usr/local/cuda/lib64 -lcudart -lcublas test.o
  </code></pre>
  
  <p>To open the MATLAB, you need to link the libstdc++ library:</p>
  
  <pre><code>LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6 matlab
  </code></pre>
  
  <p>In MATLAB, I did the following to test the code:</p>
  
  <pre><code>N = 500; x = randn(N,1);B = randn(N);A = B'*B; b = randn(N,1);    
  tic, 1/2*x'*A*x-b'*x, toc    
  tic, phi = cublas_mex_test(x,A,b),toc
  </code></pre>
  
  <p>Below is my code test.cu:</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;math.h&gt;
  #include &lt;cuda_runtime.h&gt;
  #include "cublas_v2.h"
  #include "cuda.h"
  #include "mex.h"
  #include &lt;time.h&gt;
  #define IDX2C(i,j,ld) (((j)*(ld))+(i))
  
  // this is the actual function that evaluates phi = 1/2*x'*A*x - b'*x
  int PhiEval(double *x, double *A, double *b, double &amp;phi, size_t m)
  { 
      cudaError_t cudaStat; 
      cublasStatus_t stat;
      cublasHandle_t handle;
      clock_t start, end;
  
      start = clock();
      //host data
      double *bx, *xAx;
      bx = (double*)malloc(1*sizeof(double)); 
      xAx = (double*)malloc(1*sizeof(double)); 
  
      //device data
      double* A_d, * x_d, *b_d, *Ax_d, *bx_d, *xAx_d;
      cudaStat = cudaMalloc ((void**)&amp;A_d, m*m*sizeof(double));
      cudaStat = cudaMalloc ((void**)&amp;x_d, m*sizeof(double));
      cudaStat = cudaMalloc ((void**)&amp;b_d, m*sizeof(double));
      cudaStat = cudaMalloc ((void**)&amp;Ax_d, m*sizeof(double));
      cudaStat = cudaMalloc ((void**)&amp;bx_d, 1*sizeof(double));
      cudaStat = cudaMalloc ((void**)&amp;xAx_d, 1*sizeof(double));
      if (cudaStat != cudaSuccess) {
          printf ("device memory allocation failed");
          return EXIT_FAILURE;
      }
      end = clock();
      printf ("It takes: %d clicks (%.7f seconds) to allocate the memory.\n",end-start,(double)(end-start)/CLOCKS_PER_SEC);
  
      start = clock();
      stat = cublasCreate(&amp;handle);
      if (stat != CUBLAS_STATUS_SUCCESS) {
          printf ("CUBLAS initialization failed\n");
          return EXIT_FAILURE;
      }
  
      // copy host data to device
      stat = cublasSetMatrix (m,m, sizeof(double), A, m, A_d, m);
      stat = cublasSetVector (m, sizeof(double), x, 1, x_d, 1);
      stat = cublasSetVector (m, sizeof(double), b, 1, b_d, 1);
      if (stat != CUBLAS_STATUS_SUCCESS) {
          printf ("data download failed");
          cudaFree (A_d);
          cudaFree (x_d);
          cudaFree (b_d);
          cudaFree (Ax_d);
          cudaFree(xAx_d);
          cudaFree(bx_d);
          cublasDestroy(handle);
          free(bx); free(xAx);
          return EXIT_FAILURE;
      }
      end = clock();
      printf ("It takes: %d clicks (%.7f seconds) to copy the data to GPU.\n",end-start,(double)(end-start)/CLOCKS_PER_SEC);
  
  
      start = clock();
      //calculate A*x and store the result in Ax_d
      double alpha = 1;
      double beta = 0;
      stat = cublasDgemv(handle, CUBLAS_OP_N, m,m, &amp;alpha, A_d, m, x_d, 1, &amp;beta, Ax_d, 1);
      if (stat != CUBLAS_STATUS_SUCCESS) {
          printf ("data download failed");
          cudaFree (A_d);
          cudaFree (x_d);
          cudaFree (b_d);
          cudaFree (Ax_d);
          cudaFree(xAx_d);
          cudaFree(bx_d);
          cublasDestroy(handle);
          free(bx); free(xAx);
          return EXIT_FAILURE;
      }
  
      //calculate x'*A*x and store the result in xAx_d
      stat = cublasDgemv(handle, CUBLAS_OP_T, m,1, &amp;alpha, x_d, m, Ax_d, 1, &amp;beta, xAx_d, 1);
      if (stat != CUBLAS_STATUS_SUCCESS) {
          printf ("inner product failed");
          cudaFree (A_d);
          cudaFree (x_d);
          cudaFree (b_d);
          cudaFree (Ax_d);
          cudaFree(xAx_d);
          cudaFree(bx_d);
          cublasDestroy(handle);
          free(bx); free(xAx);
          return EXIT_FAILURE;
      }
      stat = cublasGetVector (1, sizeof(double),xAx_d, 1, xAx, 1); // copy the result x'*A*x to host
  
      //calculate b'*x and store the result in bx_d
      stat = cublasDgemv(handle, CUBLAS_OP_T, m,1, &amp;alpha, b_d, m, x_d, 1, &amp;beta, bx_d, 1);
      if (stat != CUBLAS_STATUS_SUCCESS) {
          printf ("inner product failed");
          cudaFree (A_d);
          cudaFree (x_d);
          cudaFree (b_d);
          cudaFree (Ax_d);
          cudaFree(xAx_d);
          cudaFree(bx_d);
          cublasDestroy(handle);
          free(bx); free(xAx);
          return EXIT_FAILURE;
      }
      stat = cublasGetVector (1, sizeof(double),bx_d, 1, bx, 1);
      end = clock();
      printf ("It takes: %d clicks (%.7f seconds) to call functions in cublas.\n",end-start,(double)(end-start)/CLOCKS_PER_SEC);
  
      //calculate phi = 1/2*x'*A*x - b'*x
      phi = .5*xAx[0]-bx[0];
  
      start = clock();
      //free the memory
      cudaFree (A_d);
      cudaFree (x_d);
      cudaFree (b_d);
      cudaFree (Ax_d);
      cudaFree(xAx_d);
      cudaFree(bx_d);
      cublasDestroy(handle);
      free(bx); free(xAx);
  
      end = clock();
      printf ("It takes: %d clicks (%.7f seconds) to free the memory.\n",end-start,(double)(end-start)/CLOCKS_PER_SEC);
  
      return EXIT_SUCCESS;
  }
  
  /* the gateway function */
  void mexFunction( int nlhs, mxArray *plhs[],
                    int nrhs, const mxArray *prhs[])
  {
      double phi;
      double *A, *b, *x;
  
      size_t mrows,ncols;
  
      /*  check for proper number of arguments */
      if(nrhs!=3) 
          mexErrMsgIdAndTxt( "MATLAB:MinTest:invalidNumInputs",
              "Three inputs required.");
      if(nlhs!=1) 
          mexErrMsgIdAndTxt( "MATLAB:MinTest:invalidNumOutputs",
              "One output required.");
  
  
      /*  create a pointer to the input vector x */
      x = mxGetPr(prhs[0]);
  
      /*  create a pointer to the input matrix A */
      A = mxGetPr(prhs[1]);
  
      /*  create a pointer to the input vector b */
      b = mxGetPr(prhs[2]);
  
      /*  get the dimensions of the matrix input A */
      mrows = mxGetM(prhs[1]);
      ncols = mxGetN(prhs[1]);
  
      if(mrows!=ncols) 
          mexErrMsgIdAndTxt( "MATLAB:MinTest:invalidMatrixInput",
              "A has to be a square matrix");
  
      /*  call the cpp subroutine */
      PhiEval(x,A,b,phi,mrows);   
  
      plhs[0] = mxCreateDoubleScalar(phi);
  }
  </code></pre>