how to solve bugs for the matlab example "fulltosparse.c"

Question

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);

Answer 1

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 kth 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:

• At 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.
• At 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.
• When it reaches 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);
}