# HG changeset patch # User jwe # Date 1153884993 0 # Node ID 11bb9bf343a05658607fafdb763c375544519709 # Parent 6af4cea82cc7112b1944e4ca5a6cf56e92982645 [project @ 2006-07-26 03:36:33 by jwe] diff -r 6af4cea82cc7 -r 11bb9bf343a0 ChangeLog --- a/ChangeLog Tue Jul 25 19:56:00 2006 +0000 +++ b/ChangeLog Wed Jul 26 03:36:33 2006 +0000 @@ -1,3 +1,7 @@ +2006-07-25 David Bateman + + * mysparse.c: New file. + 2006-06-27 John W. Eaton * octMakefile.in (maintainer-clean distclean): Remove diff -r 6af4cea82cc7 -r 11bb9bf343a0 examples/mysparse.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/mysparse.c Wed Jul 26 03:36:33 2006 +0000 @@ -0,0 +1,116 @@ +#include "mex.h" + +void +mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + int n, m, nz; + mxArray *v; + int i; + double *pr, *pi; + double *pr2, *pi2; + int *ir, *jc; + int *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 with %d elements\n", + m, 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) = (%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])) + { + bool *pbr, *pbr2; + mexPrintf ("Matrix is %d-by-%d logical sparse matrix with %d elements\n", + m, 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 with %d elements\n", + m, 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; + } +} diff -r 6af4cea82cc7 -r 11bb9bf343a0 src/ChangeLog --- a/src/ChangeLog Tue Jul 25 19:56:00 2006 +0000 +++ b/src/ChangeLog Wed Jul 26 03:36:33 2006 +0000 @@ -1,3 +1,16 @@ +2006-07-25 David Bateman + + * mex.cc (mxArray_octave_value::get_class_id): Handle sparse. + (class mxArray_sparse): Derive from mxArray_matlab, not + mxArray_number. + (mxArray_sparse::as_octave_value): Implement function. + * ov-bool-sparse.cc (octave_sparse_bool_matrix::as_mxArray): + Implement function. + * ov-cx-sparse.cc (octave_sparse_complex_matrix::as_mxArray): + Implement function. + * ov-re-sparse.cc (octave_sparse_matrix::as_mxArray): + Implement function. + 2006-07-25 John W. Eaton * mex.cc (mxArray_struct::as_octave_value, call_mex, diff -r 6af4cea82cc7 -r 11bb9bf343a0 src/mex.cc --- a/src/mex.cc Tue Jul 25 19:56:00 2006 +0000 +++ b/src/mex.cc Wed Jul 26 03:36:33 2006 +0000 @@ -372,6 +372,13 @@ id = mxCHAR_CLASS; else if (cn == "double") id = mxDOUBLE_CLASS; + else if (cn == "sparse") + { + if (val.is_bool_type ()) + id = mxLOGICAL_CLASS; + else + id = mxDOUBLE_CLASS; + } else if (cn == "single") id = mxSINGLE_CLASS; else if (cn == "int8") @@ -1312,35 +1319,120 @@ // Matlab-style sparse arrays. -class mxArray_sparse : public mxArray_number +class mxArray_sparse : public mxArray_matlab { public: mxArray_sparse (mxClassID id_arg, int m, int n, int nzmax_arg, mxComplexity flag = mxREAL) - : mxArray_number (id_arg, m, n, flag), nzmax (nzmax_arg) + : mxArray_matlab (id_arg, m, n), nzmax (nzmax_arg) { + pr = (calloc (nzmax, get_element_size ())); + pi = (flag == mxCOMPLEX ? calloc (nzmax, get_element_size ()) : 0); ir = static_cast (calloc (nzmax, sizeof (int))); - jc = static_cast (calloc (nzmax, sizeof (int))); + jc = static_cast (calloc (n + 1, sizeof (int))); } mxArray_sparse *clone (void) const { return new mxArray_sparse (*this); } ~mxArray_sparse (void) { + mxFree (pr); + mxFree (pi); mxFree (ir); mxFree (jc); } octave_value as_octave_value (void) const { - // FIXME - abort (); - return octave_value (); + octave_value retval; + + dim_vector dv = dims_to_dim_vector (); + + switch (get_class_id ()) + { + case mxLOGICAL_CLASS: + { + bool *ppr = static_cast (pr); + + SparseBoolMatrix val (get_m (), get_n (), nzmax); + + for (int i = 0; i < nzmax; i++) + { + val.xdata(i) = ppr[i]; + val.xridx(i) = ir[i]; + } + + for (int i = 0; i < get_n () + 1; i++) + val.xcidx(i) = jc[i]; + + retval = val; + } + break; + + case mxSINGLE_CLASS: + error ("single precision data type not supported"); + break; + + case mxDOUBLE_CLASS: + { + if (pi) + { + double *ppr = static_cast (pr); + double *ppi = static_cast (pi); + + SparseComplexMatrix val (get_m (), get_n (), nzmax); + + for (int i = 0; i < nzmax; i++) + { + val.xdata(i) = Complex (ppr[i], ppi[i]); + val.xridx(i) = ir[i]; + } + + for (int i = 0; i < get_n () + 1; i++) + val.xcidx(i) = jc[i]; + + retval = val; + } + else + { + double *ppr = static_cast (pr); + + SparseMatrix val (get_m (), get_n (), nzmax); + + for (int i = 0; i < nzmax; i++) + { + val.xdata(i) = ppr[i]; + val.xridx(i) = ir[i]; + } + + for (int i = 0; i < get_n () + 1; i++) + val.xcidx(i) = jc[i]; + + retval = val; + } + } + break; + + default: + panic_impossible (); + } + + return retval; } + int is_complex (void) const { return pi != 0; } + int is_sparse (void) const { return 1; } + void *get_data (void) const { return pr; } + + void *get_imag_data (void) const { return pi; } + + void set_data (void *pr_arg) { pr = pr_arg; } + + void set_imag_data (void *pi_arg) { pi = pi_arg; } + int *get_ir (void) const { return ir; } int *get_jc (void) const { return jc; } @@ -1357,19 +1449,23 @@ int nzmax; + void *pr; + void *pi; int *ir; int *jc; mxArray_sparse (const mxArray_sparse& val) - : mxArray_number (val), nzmax (val.nzmax), + : mxArray_matlab (val), nzmax (val.nzmax), ir (static_cast (malloc (nzmax * sizeof (int)))), jc (static_cast (malloc (nzmax * sizeof (int)))) { - for (int i = 0; i < nzmax; i++) - { - ir[i] = val.ir[i]; - jc[i] = val.jc[i]; - } + int ntot = nzmax * get_element_size (); + + memcpy (pr, val.pr, ntot); + memcpy (ir, val.ir, nzmax * sizeof(int)); + memcpy (jc, val.jc, (val.get_n () + 1) * sizeof (int)); + if (pi) + memcpy (pi, val.pi, ntot); } }; diff -r 6af4cea82cc7 -r 11bb9bf343a0 src/ov-bool-sparse.cc --- a/src/ov-bool-sparse.cc Tue Jul 25 19:56:00 2006 +0000 +++ b/src/ov-bool-sparse.cc Wed Jul 26 03:36:33 2006 +0000 @@ -693,8 +693,23 @@ mxArray * octave_sparse_bool_matrix::as_mxArray (void) const { - // FIXME - return 0; + int nz = nzmax (); + mxArray *retval = new mxArray (mxLOGICAL_CLASS, rows (), columns (), + nz, mxREAL); + bool *pr = static_cast (retval->get_data ()); + int *ir = retval->get_ir (); + int *jc = retval->get_jc (); + + for (int i = 0; i < nz; i++) + { + pr[i] = matrix.data(i); + ir[i] = matrix.ridx(i); + } + + for (int i = 0; i < columns () + 1; i++) + jc[i] = matrix.cidx(i); + + return retval; } /* diff -r 6af4cea82cc7 -r 11bb9bf343a0 src/ov-cx-sparse.cc --- a/src/ov-cx-sparse.cc Tue Jul 25 19:56:00 2006 +0000 +++ b/src/ov-cx-sparse.cc Wed Jul 26 03:36:33 2006 +0000 @@ -762,8 +762,26 @@ mxArray * octave_sparse_complex_matrix::as_mxArray (void) const { - // FIXME - return 0; + int nz = nzmax (); + mxArray *retval = new mxArray (mxDOUBLE_CLASS, rows (), columns (), + nz, mxCOMPLEX); + double *pr = static_cast (retval->get_data ()); + double *pi = static_cast (retval->get_imag_data ()); + int *ir = retval->get_ir (); + int *jc = retval->get_jc (); + + for (int i = 0; i < nz; i++) + { + Complex val = matrix.data(i); + pr[i] = real (val); + pi[i] = imag (val); + ir[i] = matrix.ridx(i); + } + + for (int i = 0; i < columns() + 1; i++) + jc[i] = matrix.cidx(i); + + return retval; } /* diff -r 6af4cea82cc7 -r 11bb9bf343a0 src/ov-re-sparse.cc --- a/src/ov-re-sparse.cc Tue Jul 25 19:56:00 2006 +0000 +++ b/src/ov-re-sparse.cc Wed Jul 26 03:36:33 2006 +0000 @@ -788,8 +788,24 @@ mxArray * octave_sparse_matrix::as_mxArray (void) const { - // FIXME - return 0; + int nz = nzmax(); + int nr = rows(); + int nc = columns(); + mxArray *retval = new mxArray (mxDOUBLE_CLASS, nr, nc, nz, mxREAL); + double *pr = static_cast (retval->get_data ()); + int *ir = retval->get_ir(); + int *jc = retval->get_jc(); + + for (int i = 0; i < nz; i++) + { + pr[i] = matrix.data(i); + ir[i] = matrix.ridx(i); + } + + for (int i = 0; i < nc + 1; i++) + jc[i] = matrix.cidx(i); + + return retval; } /*