This file is in C:\Program Files\MATLAB\R2013a\extern\examples\refbook
. After mex it, I used :
aa = [1 2 3 ; 4 5 6]
fulltosparse(aa)
At the first time, the command maybe works. But try fulltosparse(aa)
for more times. You will find it will crash. Could anyone tell me why ?
mex -largeArrayDims fulltosparse.c
aa = [1 2 3; 4 5 7];
fulltosparse(aa);
fulltosparse(aa);
fulltosparse(aa);
fulltosparse(aa);
fulltosparse(aa);
fulltosparse(aa);
There appears is a bug in fulltosparse.c
with the computation of the new nzmax
when additional space for storage of non-zeros is required by the sparse matrix.
At each non-zero row i
and column j
of the m
rows and n
columns, respectively, of the full
source matrix (the k
th non-zero element), a check is made (k>=nzmax
) to ensure there is sufficient storage space in the sparse matrix data buffers (sr
, si
and irs
). If there is not enough space for more elements, the buffers are expanded via mxRealloc
with enough space for an increased nzmax
number of non-zero elements.
The problem is how it computes the new nzmax
:
mwSize oldnzmax = nzmax;
percent_sparse += 0.1;
nzmax = (mwSize)ceil((double)m*(double)n*percent_sparse);
/* make sure nzmax increases atleast by 1 */
if (oldnzmax == nzmax)
nzmax++;
The function starts with an initial percent_sparse = 0.2
, which for the 2x3
full matrix aa
corresponds to nzmax = ceil(6*0.2) = 2
, and begins looping through rows (fastest) and columns. Here's what goes wrong:
k=2
(the third element in MATLAB indexing; i=0
, j=1
), it needs to expand the buffers for the first time, and the above code runs prior to reallocation. oldnzmax
is 2. percent_sparse
is increased to 0.3, giving nzmax=ceil(6*0.3)=2
. Since oldnzmax == nzmax
, it just increments (nzmax++
) and reallocates for nzmax=3
. OK.k=3
(the fourth element; i=1
, j=1
), it goes through a similar path, increasing percent_sparse
to 0.4, computes nzmax=ceil(6*0.4)=3
and increments nzmax
to 4.k=4
(the fifth element; i=0
, j=2
), the loop iteration starts with nzmax=4
and percent_sparse=0.4
. By this point it is very clear that percent_sparse
hasn't kept up: (nzmax=4)/6 = 0.666 > percent_sparse=0.4
. Now the bug becomes obvious when percent_sparse += 0.1;
gives only 0.5 and the new nzmax=ceil(6*0.5)=3
, which is less than oldnzmax=4
!.Disaster: with k=4
and nzmax=3
, it reallocates the sparse matrix buffers (sr
, si
and irs
), and overruns both buffers:
sr[k] = pr[i]; // k=4, sr is length 3
irs[k] = i; // irs is length 3
The buffers were actually reduced in size. However, even if the test was changed to oldnzmax >= nzmax
, the buffers would still not increase in size because nzmax
is computed from a percent_sparse
that is not being incremented fast enough. There are two changes that I think are needed. First, both the test and the increment need to handle the case when nzmax<oldnzmax
:
if (oldnzmax >= nzmax)
nzmax = oldnzmax+1;
Second, just for good measure, percent_sparse
should be properly updated when the conditional is true and nzmax
gets incremented rather than computed from percent_sparse
:
if (oldnzmax >= nzmax)
{
nzmax = oldnzmax+1;
percent_sparse = (double)nzmax/((double)m*(double)n);
}