MEX Fortran Gateway and Sparse Matrices

Hello Community,

I was wondering if the following behaviour is commonplace.

Are the following two expressions the same when passed to a compiled MexFunction subroutine in fortran? Say the source code is rif.F, then

>> M = sprandsym(n, 0.5)
>> F = full(M);
>> rif(M(:,1)) % (1)
>> rif(F(:,1)) % (2)

I am obtaining different results for (1) and (2), inside MexFunction I print the copied vectors and the correct version is (2).

Illustrative example

    >> F(:,1)

    ans =

   2.088396433123533
  -0.283514990578777
                   0
   2.161892637140028
                   0
>> rif(M(:,1))
  Copied:  2.0883964331235330      -0.28351499057877672        2.1618926371400278        0.0000000000000000        0.0000000000000000
>> rif(F(:,1))
  Copied:  2.0883964331235330      -0.28351499057877672        0.0000000000000000        2.1618926371400278        0.0000000000000000

For this particular example it would seem the elements 3 and 4 were swapped, but different swaps also occur, specially with other elements in M.

Is this the expected behaviour?

Cheers!

Answer 1

[James Tursa]

Sparse matrices and full matrices are internally stored differently. For full matrices, all of the matrix values are physically stored in memory (including the 0's). For sparse matrices, only the non-zero values are stored along with row and column index information. You cannot use the same code in a mex routine to access the elements of a full and a sparse matrix. You must have different code for each. The expected result of using full code on a sparse matrix in a mex routine is either (1) wrong results, or (2) MATLAB crash.

The fix is to write specific code in your mex routine to handle sparse matrices. Options would be:

1) Write code to handle the sparse matrix directly
2) Convert the input to full before processing it (NOT recommended since this could be slow and a memory hog)
3) Throw an error if the input is sparse

Without seeing your code and knowing the variable sizes involved it is impossible to advise what approach you should take.

[Andrew]

Thanks James, this answers my query.

tl;dr

Somehow I thought M(:, 1) would expand to full before handling it to the mex subroutine.

Thanks.