changeset 5903:11bb9bf343a0

[project @ 2006-07-26 03:36:33 by jwe]
author jwe
date Wed, 26 Jul 2006 03:36:33 +0000
parents 6af4cea82cc7
children 80d3933fb8b6
files ChangeLog examples/mysparse.c src/ChangeLog src/mex.cc src/ov-bool-sparse.cc src/ov-cx-sparse.cc src/ov-re-sparse.cc
diffstat 7 files changed, 296 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2006-07-25  David Bateman  <dbateman@free.fr>
+
+	* mysparse.c: New file.
+
 2006-06-27  John W. Eaton  <jwe@octave.org>
 
 	* octMakefile.in (maintainer-clean distclean): Remove
new file mode 100644
--- /dev/null
+++ b/examples/mysparse.c
@@ -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;
+    }
+}
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,16 @@
+2006-07-25  David Bateman  <dbateman@free.fr>
+
+	* 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  <jwe@octave.org>
 
 	* mex.cc (mxArray_struct::as_octave_value, call_mex,
--- a/src/mex.cc
+++ b/src/mex.cc
@@ -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<int *> (calloc (nzmax, sizeof (int)));
-    jc = static_cast<int *> (calloc (nzmax, sizeof (int)));
+    jc = static_cast<int *> (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<bool *> (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<double *> (pr);
+	      double *ppi = static_cast<double *> (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<double *> (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<int *> (malloc (nzmax * sizeof (int)))),
       jc (static_cast<int *> (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);
   }
 };
 
--- a/src/ov-bool-sparse.cc
+++ b/src/ov-bool-sparse.cc
@@ -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<bool *> (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;
 }
 
 /*
--- a/src/ov-cx-sparse.cc
+++ b/src/ov-cx-sparse.cc
@@ -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<double *> (retval->get_data ());
+  double *pi = static_cast<double *> (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;
 }
 
 /*
--- a/src/ov-re-sparse.cc
+++ b/src/ov-re-sparse.cc
@@ -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<double *> (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;
 }
 
 /*