A.2.6 Sparse Matrices with Mex-Files

The Octave format for sparse matrices is identical to the mex format in that it is a compressed column sparse format. Also, in both implementations sparse matrices are required to be two-dimensional. The only difference of importance to the programmer is that the real and imaginary parts of the matrix are stored separately.

The mex-file interface, in addition to using mxGetM, mxGetN, mxSetM, mxSetN, mxGetPr, mxGetPi, mxSetPr, and mxSetPi, also supplies the following functions.

mwIndex *mxGetIr (const mxArray *ptr);
mwIndex *mxGetJc (const mxArray *ptr);
mwSize mxGetNzmax (const mxArray *ptr);

void mxSetIr (mxArray *ptr, mwIndex *ir);
void mxSetJc (mxArray *ptr, mwIndex *jc);
void mxSetNzmax (mxArray *ptr, mwSize nzmax);

mxGetNzmax gets the maximum number of elements that can be stored in the sparse matrix. This is not necessarily the number of nonzero elements in the sparse matrix. mxGetJc returns an array with one additional value than the number of columns in the sparse matrix. The difference between consecutive values of the array returned by mxGetJc define the number of nonzero elements in each column of the sparse matrix. Therefore,

mwSize nz, n;
mwIndex *Jc;
mxArray *m;
…
n = mxGetN (m);
Jc = mxGetJc (m);
nz = Jc[n];

returns the actual number of nonzero elements stored in the matrix in nz. As the arrays returned by mxGetPr and mxGetPi only contain the nonzero values of the matrix, we also need a pointer to the rows of the nonzero elements, and this is given by mxGetIr. A complete example of the use of sparse matrices in mex-files is given by the file mysparse.c shown below.

#include "mex.h"

void
mexFunction (int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
{
  mwSize m, n, nz;
  mxArray *v;
  mwIndex i;
  double *pr, *pi;
  double *pr2, *pi2;
  mwIndex *ir, *jc;
  mwIndex *ir2, *jc2;

  if (nrhs != 1 || ! mxIsSparse (prhs[0]))
    mexErrMsgTxt ("ARG1 must be a sparse matrix");

  m = mxGetM (prhs[0]);
  n = mxGetN (prhs[0]);
  nz = mxGetNzmax (prhs[0]);

  if (mxIsComplex (prhs[0]))
    {
      mexPrintf ("Matrix is %d-by-%d complex sparse matrix", m, n);
      mexPrintf (" with %d elements\n", nz);

      pr = mxGetPr (prhs[0]);
      pi = mxGetPi (prhs[0]);
      ir = mxGetIr (prhs[0]);
      jc = mxGetJc (prhs[0]);

      i = n;
      while (jc[i] == jc[i-1] && i != 0) i--;

      mexPrintf ("last nonzero element (%d, %d) = (%g, %g)\n",
                 ir[nz-1]+ 1, i, pr[nz-1], pi[nz-1]);

      v = mxCreateSparse (m, n, nz, mxCOMPLEX);
      pr2 = mxGetPr (v);
      pi2 = mxGetPi (v);
      ir2 = mxGetIr (v);
      jc2 = mxGetJc (v);

      for (i = 0; i < nz; i++)
        {
          pr2[i] = 2 * pr[i];
          pi2[i] = 2 * pi[i];
          ir2[i] = ir[i];
        }
      for (i = 0; i < n + 1; i++)
        jc2[i] = jc[i];

      if (nlhs > 0)
        plhs[0] = v;
    }
  else if (mxIsLogical (prhs[0]))
    {
      mxLogical *pbr, *pbr2;
      mexPrintf ("Matrix is %d-by-%d logical sparse matrix", m, n);
      mexPrintf (" with %d elements\n", nz);

      pbr = mxGetLogicals (prhs[0]);
      ir = mxGetIr (prhs[0]);
      jc = mxGetJc (prhs[0]);

      i = n;
      while (jc[i] == jc[i-1] && i != 0) i--;
      mexPrintf ("last nonzero element (%d, %d) = %d\n",
                 ir[nz-1]+ 1, i, pbr[nz-1]);

      v = mxCreateSparseLogicalMatrix (m, n, nz);
      pbr2 = mxGetLogicals (v);
      ir2 = mxGetIr (v);
      jc2 = mxGetJc (v);

      for (i = 0; i < nz; i++)
        {
          pbr2[i] = pbr[i];
          ir2[i] = ir[i];
        }
      for (i = 0; i < n + 1; i++)
        jc2[i] = jc[i];

      if (nlhs > 0)
        plhs[0] = v;
    }
  else
    {
      mexPrintf ("Matrix is %d-by-%d real sparse matrix", m, n);
      mexPrintf (" with %d elements\n", nz);

      pr = mxGetPr (prhs[0]);
      ir = mxGetIr (prhs[0]);
      jc = mxGetJc (prhs[0]);

      i = n;
      while (jc[i] == jc[i-1] && i != 0) i--;
      mexPrintf ("last nonzero element (%d, %d) = %g\n",
                 ir[nz-1]+ 1, i, pr[nz-1]);

      v = mxCreateSparse (m, n, nz, mxREAL);
      pr2 = mxGetPr (v);
      ir2 = mxGetIr (v);
      jc2 = mxGetJc (v);

      for (i = 0; i < nz; i++)
        {
          pr2[i] = 2 * pr[i];
          ir2[i] = ir[i];
        }
      for (i = 0; i < n + 1; i++)
        jc2[i] = jc[i];

      if (nlhs > 0)
        plhs[0] = v;
    }
}

A sample usage of mysparse is

sm = sparse ([1, 0; 0, pi]);
mysparse (sm)
⇒
Matrix is 2-by-2 real sparse matrix with 2 elements
last nonzero element (2, 2) = 3.14159

© 1996–2020 John W. Eaton
Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies.
Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one.
Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions.
https://octave.org/doc/v5.2.0/Sparse-Matrices-with-Mex_002dFiles.html