Info Catalog octave: Structures with Mex-Files octave: Mex-Files octave: Calling Other Functions in Mex-Files

octave: Sparse Matrices with Mex-Files

 
 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
 

Info Catalog octave: Structures with Mex-Files octave: Mex-Files octave: Calling Other Functions in Mex-Files