Mercurial > octave-dspies
view examples/mysparse.c @ 19006:2e0613dadfee draft
All calls to "find" use the same generic implementation (bug #42408, 42421)
* find.cc: Rewrite.
Move generic "find" logic to find.h
(Ffind) : Changed calls to find_nonzero_elem_idx to find_templated
Added unit test for bug #42421
* Array.cc (and .h) (Array::find): Deleted function. Replaced with find::find(Array)
from find.h
* Array.h: Added typedef for array_iterator (in nz-iterators.h) as
Array::iter_type
* DiagArray2.h: Added typedef for diag_iterator (in nz-iterators.h) as
DiagArray2::iter_type
* PermMatrix.h: Added typedef for perm_iterator (in nz-iterators.h) as
PermMatrix::iter_type
Also added typedef for bool as PermMatrix::element_type
(not octave_idx_type)
Added an nnz() function (which is an alias for perm_length) and a
perm_elem(i) function for retrieving the ith element of the permutation
* Sparse.h: Added typedef for sparse_iterator (in nz-iterators.h) as
Sparse::iter_type
Added a short comment documenting the the argument to the numel
function
* idx-vector.cc (idx_vector::idx_mask_rep::as_array): Changed Array.find to
find::find(Array) (in find.h)
* (new file) find.h
* (new file) interp-idx.h: Simple methods for converting between interpreter
index type and internal octave_idx_type/row-col pair
* (new file) min-with-nnz.h: Fast methods for taking an arbitrary matrix M and
an octave_idx_type n and finding min(M.nnz(), n)
* (new file) nz-iterators.h: Iterators for traversing (in column-major order)
the nonzero elements of any array or matrix backwards or forwards
* (new file) direction.h: Generic methods for simplifying code has to deal with
a "backwards or forwards" template argument
* build-sparse-tests.sh: Removed 5-return-value calls to "find" in unit-tests;
Admittedly this commit breaks this "feature" which was undocumented and only
partially supported to begin with (ie never worked for full matrices,
permutation matrices, or diagonal matrices)
author | David Spies <dnspies@gmail.com> |
---|---|
date | Tue, 17 Jun 2014 16:41:11 -0600 |
parents | 9ac2357f19bc |
children |
line wrap: on
line source
#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; } }