<html lang="en">
<head>
<title>Sparse Matrices with Mex-Files - GNU Octave</title>
<meta http-equiv="Content-Type" content="text/html">
<meta name="description" content="GNU Octave">
<meta name="generator" content="makeinfo 4.13">
<link title="Top" rel="start" href="index.html#Top">
<link rel="up" href="Mex_002dFiles.html#Mex_002dFiles" title="Mex-Files">
<link rel="prev" href="Structures-with-Mex_002dFiles.html#Structures-with-Mex_002dFiles" title="Structures with Mex-Files">
<link rel="next" href="Calling-Other-Functions-in-Mex_002dFiles.html#Calling-Other-Functions-in-Mex_002dFiles" title="Calling Other Functions in Mex-Files">
<link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
<meta http-equiv="Content-Style-Type" content="text/css">
<style type="text/css"><!--
pre.display { font-family:inherit }
pre.format { font-family:inherit }
pre.smalldisplay { font-family:inherit; font-size:smaller }
pre.smallformat { font-family:inherit; font-size:smaller }
pre.smallexample { font-size:smaller }
pre.smalllisp { font-size:smaller }
span.sc { font-variant:small-caps }
span.roman { font-family:serif; font-weight:normal; }
span.sansserif { font-family:sans-serif; font-weight:normal; }
--></style>
</head>
<body>
<div class="node">
<a name="Sparse-Matrices-with-Mex-Files"></a>
<a name="Sparse-Matrices-with-Mex_002dFiles"></a>
<p>
Next: <a rel="next" accesskey="n" href="Calling-Other-Functions-in-Mex_002dFiles.html#Calling-Other-Functions-in-Mex_002dFiles">Calling Other Functions in Mex-Files</a>,
Previous: <a rel="previous" accesskey="p" href="Structures-with-Mex_002dFiles.html#Structures-with-Mex_002dFiles">Structures with Mex-Files</a>,
Up: <a rel="up" accesskey="u" href="Mex_002dFiles.html#Mex_002dFiles">Mex-Files</a>
<hr>
</div>
<h4 class="subsection">A.2.6 Sparse Matrices with Mex-Files</h4>
<p>The Octave format for sparse matrices is identical to the mex format in
that it is a compressed column sparse format. Also in both, sparse
matrices are required to be two-dimensional. The only difference is that
the real and imaginary parts of the matrix are stored separately.
<p>The mex-file interface, as well as using <code>mxGetM</code>, <code>mxGetN</code>,
<code>mxSetM</code>, <code>mxSetN</code>, <code>mxGetPr</code>, <code>mxGetPi</code>,
<code>mxSetPr</code> and <code>mxSetPi</code>, the mex-file interface supplies the
functions
<pre class="example"> 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);
</pre>
<p class="noindent"><code>mxGetNzmax</code> gets the maximum number of elements that can be stored
in the sparse matrix. This is not necessarily the number of non-zero
elements in the sparse matrix. <code>mxGetJc</code> 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
<code>mxGetJc</code> define the number of non-zero elements in each column of
the sparse matrix. Therefore
<pre class="example"> mwSize nz, n;
mwIndex *Jc;
mxArray *m;
...
n = mxGetN (m);
Jc = mxGetJc (m);
nz = Jc[n];
</pre>
<p class="noindent">returns the actual number of non-zero elements stored in the matrix in
<code>nz</code>. As the arrays returned by <code>mxGetPr</code> and <code>mxGetPi</code>
only contain the non-zero values of the matrix, we also need a pointer
to the rows of the non-zero elements, and this is given by
<code>mxGetIr</code>. A complete example of the use of sparse matrices in
mex-files is given by the file <samp><span class="file">mysparse.c</span></samp> as seen below
<pre class="example"><pre class="verbatim"> #include "mex.h"
void
mexFunction (int nlhs, mxArray *plhs[], int nrhs,
const mxArray *prhs[])
{
mwSize n, m, nz;
mxArray *v;
mwIndex i;
double *pr, *pi;
double *pr2, *pi2;
mwIndex *ir, *jc;
mwIndex *ir2, *jc2;
if (nrhs != 1 || ! mxIsSparse (prhs[0]))
mexErrMsgTxt ("expects 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 non-zero element (%d, %d) =",
ir[nz-1]+ 1, i);
mexPrintf (" (%g, %g)\n", 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]))
{
bool *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 non-zero 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 non-zero 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;
}
}
</pre>
</pre>
</body></html>
