MATLAB Answers

Sending big sparse matrix to C through mex: matlab crashing

18 views (last 30 days)
Dushyant
Dushyant on 9 Feb 2015
Edited: Dushyant on 10 Feb 2015
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);
}

  0 Comments

Sign in to comment.

Accepted Answer

James Tursa
James Tursa on 9 Feb 2015
Edited: James Tursa on 9 Feb 2015
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).

  3 Comments

Dushyant
Dushyant on 9 Feb 2015
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
James Tursa on 9 Feb 2015
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
Dushyant on 10 Feb 2015
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:
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?

Sign in to comment.

More Answers (0)

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!