Category Archives: Matrix

Implementation mistake causes memory leakage when build Sparse Matrix in mex C/C++ to interact with MATLAB

Let’s first take a look at 6 Lessons From Dropbox – One Million Files Saved Every 15 Minutes, we know, when the product scales, engineers tend to use more computational efficient (C/C++) languages to compensate the side effect of rapid prototyping/development (Python) language.

Similar thing happens in academic prototyping, MATLAB is convenient, and more worthwhile to be used to explore novel idea and save time. However, some important parts are required to implemented in more efficient way to make the idea more solid. Therefore, programming and interacting programming languages is a handy skill set to computer science students.  When we deal with those situation, memory management and transferring is always pain.

In this article, I program in C with MEX to talk to MATLAB. I mistakenly transfer the sparse matrix format with redundant allocations so it blows up the memory (PS: I deal with large scale computing problem).

Do not use this memory Blow-Up version code:

plhs[0] = mxCreateSparse( c_stamp->M, c_stamp->N, c_stamp->NNZ, mxREAL);
cPr = (double *)mxMalloc(sizeof(double) * (c_stamp->NNZ));
cIr = (mwIndex *)mxMalloc(sizeof(mwIndex) * (c_stamp->NNZ));
cJc = (mwIndex *)mxMalloc(sizeof(mwIndex) * (c_stamp->N + 1));

//copy memory
memcpy( cPr, c_stamp->Pr, sizeof(double) * (c_stamp->NNZ));
memcpy( cIr, (mwIndex *)(c_stamp->Ir), sizeof(mwIndex) * (c_stamp->NNZ));
memcpy( cJc, (mwIndex *)(c_stamp->Jc), sizeof(mwIndex) * (c_stamp->N + 1));

mxSetPr(plhs[0], cPr);
mxSetIr(plhs[0], cIr);
mxSetJc(plhs[0], cJc);

Instead, this version is right

plhs[0] = mxCreateSparse( c_stamp->M, c_stamp->N, c_stamp->NNZ, mxREAL);
cPr = (double *) mxGetPr(plhs[0]);
memcpy(cPr, c_stamp->Pr, sizeof(double)*(c_stamp->NNZ));
cIr = (mwIndex *) mxGetIr(plhs[0]);
memcpy( cIr, (mwIndex *)(c_stamp->Ir), sizeof(mwIndex) * (c_stamp->NNZ));
cJc = (mwIndex *) mxGetJc(plhs[0]);
memcpy( cJc, (mwIndex *)(c_stamp->Jc), sizeof(mwIndex) * (c_stamp->N + 1));

where cPr, cIr, cJc are defined before. c_stamp is a sparse matrix I want to wrap from C/C++ to MATLAB via mex.


Julia @ Ubuntu 12.04

Julia has some issues when I compiled in Ubuntu 12.04. I modify a little on top of

1.)  at terminals

sudo add-apt-repository ppa:ubuntu-toolchain-r/test

2.) Then install gcc 4.8,  g++ 4.8, gfortran-4.8

sudo apt-get update; sudo apt-get install gcc-4.8 g++-4.8 gfortran-4.8

3.) Once installed, run following commands one by one to use gcc 4.8 instead of previous version.

sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-4.8 20
sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-4.8 20
sudo update-alternatives --install /usr/bin/gfortran gfortran /usr/bin/gfortran-4.8 20

Besides stiffness problem, the explicit methods are still not better than implicit methods?

Considering the simple dynamical system as follows

{\bf C}{\bf \dot{x}}(t)={\bf G}{\bf }x(t)+{\bf B}{\bf u}(t)

where \bf{C} and \bf{G} are sparse matrices, {\bf u}(t) is a input and \bf{B} is incident matrix which poses {\bf u}(t) in the system.

In my thought, the explicit methods work if there is not \bf{C} on the left hand side,

{\bf \dot{x}}(t)={\bf G}{\bf x}(t)+{\bf B}{\bf u}(t)

e.g. solve it by forward Euler method and ignore the input for simplicity (\bf{\dot{x}}(t)={\bf G}{\bf x}(t))), then

{\bf x}(t+h)={\bf x}(t)+h{\bf G}{\bf x}(t), where h is the discretized time step.

However, when this is not true, the explicit method still need to solve linear system, because

{\bf C}{\bf x}(t+h)={\bf C}{\bf x}(t)+h{\bf G}{\bf x}(t) or say it requires inversion of matrix {\bf x}(t+h)={\bf x}(t)+h{\bf C}^{-1}{\bf G}{\bf x}(t)

Any comments?