Sending big sparse matrix to C through mex: matlab crashing

I have written a mex code which sends a scalar and a matrix to C-code from matlab code. It works fine with smaller matrix. However, when I try to pass big sparse matrix (size ~ 8448 x 3264), matlab crashes with following error:

I get the following error: Matlab has encountered an internal problem and needs to close.

                      *Detailed error report*
  ------------------------------------------------------------------------
         Segmentation violation detected at Mon Feb  9 13:21:48 2015
  ------------------------------------------------------------------------

Configuration:
  Crash Decoding     : Disabled
  Current Visual     : None
  Default Encoding   : UTF-8
  GNU C Library      : 2.19 stable
  MATLAB Architecture: glnxa64
  MATLAB Root        : /usr/local/MATLAB/R2014b
  MATLAB Version     : 8.4.0.150421 (R2014b)
  Operating System   : Linux 3.13.0-37-generic #64-Ubuntu SMP Mon Sep 22 21:28:38 UTC 2014 x86_64
  Processor ID       : x86 Family 6 Model 44 Stepping 2, GenuineIntel
  Software OpenGL    : 1
  Virtual Machine    : Java 1.7.0_11-b21 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
  Window System      : No active display

Fault Count: 1

Abnormal termination:
Segmentation violation

Register State (from fault):
  RAX = 00007f011f000000  RBX = 0000000000000001
  RCX = 0000000000260fe0  RDX = 00007f0162197000
  RSP = 00007f024fffb4f0  RBP = 00007f024fffb4f0
  RSI = 00007f011ed9f020  RDI = 00007f0161f36020

     R8 = 00007f011ed9f010   R9 = 0000000000000000
    R10 = 0000000000000022  R11 = 0000000048000001
    R12 = 00007f024fffbaf0  R13 = 00007f01c63b57f0
    R14 = 00007f024fffbb00  R15 = 00007f024fffbb00

    RIP = 00007f01618d491e  EFL = 0000000000010206

     CS = 0033   FS = 0000   GS = 0000

Stack Trace (from fault):
[  0] 0x00007f01618d491e /home/dkumar/Mex_Codes_DKU/Matlab_Calling_C/Test_2/mexcallingmatlab.mexa64+00002334
[  1] 0x00007f01618d4a8e /home/dkumar/Mex_Codes_DKU/Matlab_Calling_C/Test_2/mexcallingmatlab.mexa64+00002702 mexFunction+00000325
.......

Update: Thanks to Mr. Tursa, I could sort out my problem. I am just updating this section with working code

Here is my matlab code:

% Create system matrix (size 8448 x 3264) smA_System = ConstructSystemMatrix();

tic

    x = 9;
    y = ones(3);

    % This works fine
    z = mexcallingmatlab(x, y);

    % This gives error
    z = mexcallingmatlab(x, smA_System);

toc

UPDATED PART:

Here is my mex code that works

#include "mex.h"
#include "matrix.h"   // definition of Matlab sparse matrix

#include<stdio.h>
#include "/home/dkumar/libtsnnls-2.3.3/tsnnls/taucs_basic/taucs.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  double *y, *z, x;
  
  /* Declare variable */
    mwSize mrows,ncols;
  int status;
    mwSize nzmax;
    mwIndex *irs,*jcs,j,k;

  if (nrhs != 2) 
    mexErrMsgTxt("Two inputs required.");
  if (nlhs != 1) 
    mexErrMsgTxt("One output required.");

  /* Check to make sure the first input argument is a scalar. */
  if (!mxIsDouble(prhs[0]) ||   mxIsComplex(prhs[0]) ||  mxGetN(prhs[0])*mxGetM(prhs[0]) != 1) {
    mexErrMsgTxt("Input x must be a scalar.");
  }

  /* Get the scalar input x. */
  x = mxGetScalar(prhs[0]);

  /* Create a pointer to the input matrix y. */
  y = mxGetPr(prhs[1]);

       /* Check data type dimensions of the matrix input y. */
      if (!(mxIsDouble(prhs[1]))){
          mexErrMsgIdAndTxt( "MATLAB:DKU", "Input argument must be of type double.");
      }

      if (mxGetNumberOfDimensions(prhs[1]) != 2){
          mexErrMsgIdAndTxt( "MATLAB:DKU", "Input argument must be two dimensional\n");
      }

   /* Get the size and pointers to input data */
    mrows =mxGetM(prhs[1]);
    ncols = mxGetN(prhs[1]);

  // Verify the matrix y is infact sparse
  if (!(mxIsSparse(prhs[1]))){
        mexErrMsgIdAndTxt( "MATLAB: DKU", "Input matrix is not sparse\n");
  }else{    
    sr  = mxGetPr(prhs[1]);
    jcs = mxGetJc(prhs[1]);    
    irs = mxGetIr(prhs[1]);
  }
  
   // Now converting sparse matrix to TAUCS format
     taucs_ccs_matrix *A;
    A = (taucs_ccs_matrix*)malloc(sizeof(taucs_ccs_matrix));
  A->n = ncols;
  A->m = mrows;  
  A->flags = TAUCS_DOUBLE;

      // Allocating spaces
    int nnz = jcs[ncols]; 
    A->colptr     = (int*)mxMalloc(sizeof(int)*(A->n+1));
    A->rowind   = (int*)mxMalloc(sizeof(int)*nnz);
    A->values.d   = (double*)mxMalloc(sizeof(taucs_double)*nnz);  
    
    int icolOffset = 0;  // Matlab "SPARSE" indices start with 0 vs 0 based in C  
    
    A->colptr = jcs-icolOffset;
    A->rowind = irs-icolOffset;
    A->values.d = sr;

  /* Create a C pointer to a copy of the output matrix. */
  memcpy(mxGetPr(plhs[0]), &A, sizeof(A)); 

  
  //Freeing spaces
    mxFree(A->colptr);
    mxFree(A->rowind);
    mxFree(A->values.d);  
  
}

Answer 1

[James Tursa]

Sparse matrices are not stored in memory the same way that full matrices are. In particular, the 0's are not physically stored. E.g., in a 100 x 100 full matrix there are 10,000 physical double values stored in memory. But in a 100 x 100 sparse matrix there could only be a fraction of those values physically stored in memory. If you pass a sparse matrix into a mex routine and then attempt to loop over all of the values as if they were all physically present in memory, you will likely run off the end of valid memory and get a seg fault (as you have discovered).

The fix: You need to have separate code to deal with the special sparse matrix internal format (look at some examples in the doc) ... you can't use your double loops as shown in your code. And if you want a sparse matrix output you will need to to create it with the mxCreateSparse function (not mxCreateDoubleMatrix). E.g., here is code that will tell you how many elements of the sparse matrix are physically stored in memory:

mwIndex numel;
mwIndex *Jc;
size_t N;
N = mxGetN(prhs[1]); % Number of columns
Jc = nxGetJc(prhs[1]); % Number of elements physically stored in each column (pointer to)
numel = Jc[N]; // Number of elements physically stored in memory (last value in Jc array)

The numel value could be any number from 0 to mxGetM(prhs[1]) * mxGetN(prhs[1]).

Bottom line is you need to either disallow sparse matrices as inputs to your mex routine (check for them with mxIsSparse function), or supply completely new code to deal with them (e.g. review the mex sparse code examples in the doc).

[Dushyant]

Mr. Tursa, Thanks a lot.

Another request: I have updated my question with my new "mex" code. This code passes a scalar and a sparse matrix from Matlab and convert sparse matrix into another format TAUCS (C-lib) which is also compressed column form (more info on page 12 of info about TAUCS). TAUCS format is also a column major format and appears identical to Matlab sparse format, except some flag etc.

Conversion part is in the end of the code.

Could you please check if I am doing some mistake?

PS: Of course, I am just asking your opion and I would not held you accountable if something goes wrong :).

[James Tursa]

Here is the struct definition from the pdf file in the link:

typedef struct {
int n; // number of columns
int m; // number of rows
int flags; // see below
int* colptr; // pointers to where columns begin in rowind and values
             // 0-based, length is (n+1)
int* rowind; // row indices, 0-based
union {
void* v;
taucs_double* d;
taucs_single* s;
taucs_dcomplex* z;
taucs_scomplex* c; }
values; numerical values
} taucs_ccs_matrix;

The code you list above, repeated here for convenience, is:

  int nnz = jcs[ncols]; 
  A->colptr   = (int*)malloc(sizeof(int)*(A->n+1));
  A->rowind   = (int*)malloc(sizeof(int)*nnz);
  A->values.d   = (double*)malloc(sizeof(taucs_double)*nnz);

  int icolOffset = 1;//Matlab indices start with 1 vs 0 based in C  
  
  A->colptr = jcs-icolOffset;
  A->rowind = irs-icolOffset;
  A->values.d = sr;

Unfortunately, there is not enough information in the pdf file to tell me what exactly should be in the colptr array. In MATLAB, the first element of Jc is always 0, and the remaining values in Jc are the accumulated number of elements in the array up to and including the column. E.g., Jc[1] is the total number of elements in 1st column. Jc[2] is the total number of elements in 1st & 2nd columns. Jc[3] is the total number of elements in 1st & 2nd & 3rd columns. Etc. You can think of these values as "pointers" into the Ir and Pr arrays for where each column starts (same as TAUCS???). Also, the row indices are 0-based. But in the TAUCS doc there is no explicit mention (that I can find) of what should be in colptr. From the text, I would guess that it is the same as MATLAB but I can't be sure.

Regardless, your code is definitely wrong. In particular the icolOffset needs to go (indices in the MATLAB sparse format are 0-based). Assuming mwIndex is basically an int or unsigned int (are you running 32-bit MATLAB?), you might be able to get away with simple assignments. Also, your malloc's make no sense and in fact generate memory leaks (you allocate memory and then wipe out the pointers to this memory downstream). E.g.,

  // A->colptr   = (int*)malloc(sizeof(int)*(A->n+1));
  // A->rowind   = (int*)malloc(sizeof(int)*nnz);
  // A->values.d   = (double*)malloc(sizeof(taucs_double)*nnz);

  // int icolOffset = 1;

  A->colptr = jcs;
  A->rowind = irs;
  A->values.d = sr;

I would give it a try and see what happens. (Unless you are running 64-bit MATLAB, in which case you will need the malloc's and copies of Ir and Jc).

Finally, since this is a mex routine, I would advise using mxMalloc instead of malloc for all of your memory allocations.

[Dushyant]

Thanks, Mr. Tursa.

I am using 64 bit Matlab R2014b.

Even lack of proper documentation of TAUCS is making it tougher for me. Anyway, after extensive search online, I found a C++ code which converts row-major C++ format to column major TAUCS format:

Conversion code in CPP

Please have a look at code segment starting at following comment

// 2b) converting SparseMatrix to the TAUCS CCS

Also, I modified my mex-code as you instructed.

Since I am not that good in C, I would appreciate if you could have a look at my mex code again and compare it with other C++ code.

Last part of my code: Do I need to assign memory to plhs[0] before coping TAUCS matrix into it?