Mercurial > hg > octave-max

changeset 3864:e78705239df5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2001-11-16 12:56:20 by jwe]
author jwe
date Fri, 16 Nov 2001 12:56:22 +0000
parents 55c8eee153cb
children 4592f72b7317
files liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc liboctave/mx-inlines.cc src/ChangeLog src/data.cc
diffstat 6 files changed, 149 insertions(+), 348 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2444,195 +2444,42 @@
 ComplexMatrix
 ComplexMatrix::cumprod (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  ComplexMatrix retval (nr, nc);
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  for (int i = 0; i < nr; i++)
-	    {
-	      Complex prod = elem (i, 0);
-	      for (int j = 0; j < nc; j++)
-		{
-		  retval.elem (i, j) = prod;
-		  if (j < nc - 1)
-		    prod *= elem (i, j+1);
-		}
-	    }
-	}
-      else
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      Complex prod = elem (0, j);
-	      for (int i = 0; i < nr; i++)
-		{
-		  retval.elem (i, j) = prod;
-		  if (i < nr - 1)
-		    prod *= elem (i+1, j);
-		}
-	    }
-	}
-    }
-  return retval;
+  MX_CUMMULATIVE_OP (ComplexMatrix, Complex, *=);
 }
 
 ComplexMatrix
 ComplexMatrix::cumsum (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  ComplexMatrix retval (nr, nc);
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  for (int i = 0; i < nr; i++)
-	    {
-	      Complex sum = elem (i, 0);
-	      for (int j = 0; j < nc; j++)
-		{
-		  retval.elem (i, j) = sum;
-		  if (j < nc - 1)
-		    sum += elem (i, j+1);
-		}
-	    }
-	}
-      else
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      Complex sum = elem (0, j);
-	      for (int i = 0; i < nr; i++)
-		{
-		  retval.elem (i, j) = sum;
-		  if (i < nr - 1)
-		    sum += elem (i+1, j);
-		}
-	    }
-	}
-    }
-  return retval;
+  MX_CUMMULATIVE_OP (ComplexMatrix, Complex, +=);
 }
 
 ComplexMatrix
 ComplexMatrix::prod (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  ComplexMatrix retval;
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize(nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 1.0;
-	      for (int j = 0; j < nc; j++)
-		retval.elem (i, 0) *= elem (i, j);
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 1.0;
-	      for (int i = 0; i < nr; i++)
-		retval.elem (0, j) *= elem (i, j);
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1,1);
-      retval.elem (0,0) = 1.0;
-    }
-  return retval;
+  MX_REDUCTION_OP (ComplexMatrix, *=, 1.0, 1.0);
 }
 
 ComplexMatrix
 ComplexMatrix::sum (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  ComplexMatrix retval;
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize (nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 0.0;
-	      for (int j = 0; j < nc; j++)
-		retval.elem (i, 0) += elem (i, j);
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 0.0;
-	      for (int i = 0; i < nr; i++)
-		retval.elem (0, j) += elem (i, j);
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-    }
-  return retval;
+  MX_REDUCTION_OP (ComplexMatrix, +=, 0.0, 0.0);
 }
 
 ComplexMatrix
 ComplexMatrix::sumsq (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  ComplexMatrix retval;
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize (nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 0.0;
-	      for (int j = 0; j < nc; j++)
-		{
-		  Complex d = elem (i, j);
-		  retval.elem (i, 0) += d * conj (d);
-		}
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 0.0;
-	      for (int i = 0; i < nr; i++)
-		{
-		  Complex d = elem (i, j);
-		  retval.elem (0, j) += d * conj (d);
-		}
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-    }
-
-  return retval;
+#define ROW_EXPR \
+  Complex d = elem (i, j); \
+  retval.elem (i, 0) += d * conj (d)
+
+#define COL_EXPR \
+  Complex d = elem (i, j); \
+  retval.elem (0, j) += d * conj (d)
+
+  MX_BASE_REDUCTION_OP (ComplexMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0);
+
+#undef ROW_EXPR
+#undef COL_EXPR
 }
 
 ComplexColumnVector
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,19 @@
+2001-11-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* mx-inlines.cc (MX_CUMMULATIVE_OP): New macro.
+	* CMatrix.cc (ComplexMatrix::cumprod, ComplexMatrix::cumsum): Use it.
+	* dMatrix.cc (Matrix::cumprod, Matrix::cumsum): Likewise.
+
+	* mx-inlines.cc (MX_REDUCTION_OP, MX_REDUCTION_OP_COL_EXPR,
+	MX_REDUCTION_OP_ROW_EXPR): New macros.
+	* dMatrix.cc (Matrix::prod, Matrix::sum): Use MX_REDUCTION_OP.
+	* CMatrix.cc (ComplexMatrix::prod, ComplexMatrix::sum): Likewise.
+
+	* mx-inlines.cc (MX_BASE_REDUCTION_OP): New macro.
+	DIM == -1 now means no orientation for vector sums.
+	* dMatrix.cc (ComplexMatrix::sumsq): Use it.
+	* CMatrix.cc (ComplexMatrix::sumsq): Likewise.
+
 2001-11-08  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Range.cc (Range::nelem_internal): Special case ranges that must
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -1991,203 +1991,42 @@
 Matrix
 Matrix::cumprod (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  Matrix retval (nr, nc);
-
-  if (nr > 0 && nc >0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  for (int i = 0; i < nr; i++)
-	    {
-	      double prod = elem (i, 0);
-	      for (int j = 0; j < nc; j++)
-		{
-		  retval.elem (i, j) = prod;
-		  if (j < nc - 1)
-		    prod *= elem (i, j+1);
-		}
-	    }
-	}
-      else
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      double prod = elem (0, j);
-	      for (int i = 0; i < nr; i++)
-		{
-		  retval.elem (i, j) = prod;
-		  if (i < nr - 1)
-		    prod *= elem (i+1, j);
-		}
-	    }
-	}
-    }
-  return retval;
+  MX_CUMMULATIVE_OP (Matrix, double, *=);
 }
 
 Matrix
 Matrix::cumsum (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  Matrix retval (nr, nc);
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  for (int i = 0; i < nr; i++)
-	    {
-	      double sum = elem (i, 0);
-	      for (int j = 0; j < nc; j++)
-		{
-		  retval.elem (i, j) = sum;
-		  if (j < nc - 1)
-		    sum += elem (i, j+1);
-		}
-	    }
-	}
-      else
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      double sum = elem (0, j);
-	      for (int i = 0; i < nr; i++)
-		{
-		  retval.elem (i, j) = sum;
-		  if (i < nr - 1)
-		    sum += elem (i+1, j);
-		}
-	    }
-	}
-    }
-  return retval;
+  MX_CUMMULATIVE_OP (Matrix, double, +=);
 }
 
 Matrix
 Matrix::prod (int dim) const
 {
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize(nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 1.0;
-	      for (int j = 0; j < nc; j++)
-		retval.elem (i, 0) *= elem (i, j);
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 1.0;
-	      for (int i = 0; i < nr; i++)
-		retval.elem (0, j) *= elem (i, j);
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 1.0;
-    }
-
-  return retval;
+  MX_REDUCTION_OP (Matrix, *=, 1.0, 1.0);
 }
 
 Matrix
 Matrix::sum (int dim) const
 {
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize (nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 0.0;
-	      for (int j = 0; j < nc; j++)
-		retval.elem (i, 0) += elem (i, j);
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 0.0;
-	      for (int i = 0; i < nr; i++)
-		retval.elem (0, j) += elem (i, j);
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-    }
-  return retval;
+  MX_REDUCTION_OP (Matrix, +=, 0.0, 0.0);
 }
 
 Matrix
 Matrix::sumsq (int dim) const
 {
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize (nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 0.0;
-	      for (int j = 0; j < nc; j++)
-		{
-		  double d = elem (i, j);
-		  retval.elem (i, 0) += d * d;
-		}
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 0.0;
-	      for (int i = 0; i < nr; i++)
-		{
-		  double d = elem (i, j);
-		  retval.elem (0, j) += d * d;
-		}
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-    }
-  return retval;
+#define ROW_EXPR \
+  double d = elem (i, j); \
+  retval.elem (i, 0) += d * d
+
+#define COL_EXPR \
+  double d = elem (i, j); \
+  retval.elem (0, j) += d * d
+
+  MX_BASE_REDUCTION_OP (Matrix, ROW_EXPR, COL_EXPR, 0.0, 0.0);
+
+#undef ROW_EXPR
+#undef COL_EXPR
 }
 
 Matrix
--- a/liboctave/mx-inlines.cc
+++ b/liboctave/mx-inlines.cc
@@ -227,6 +227,101 @@
 OP_DUP_FCN (imag, mx_inline_imag_dup, double,  Complex)
 OP_DUP_FCN (conj, mx_inline_conj_dup, Complex, Complex)
 
+// Avoid some code duplication.  Maybe we should use templates.
+
+#define MX_CUMMULATIVE_OP(RET_TYPE, ELT_TYPE, OP) \
+ \
+  int nr = rows (); \
+  int nc = cols (); \
+ \
+  RET_TYPE retval (nr, nc); \
+ \
+  if (nr > 0 && nc > 0) \
+    { \
+      if ((nr == 1 && dim == -1) || dim == 1) \
+	{ \
+	  for (int i = 0; i < nr; i++) \
+	    { \
+	      ELT_TYPE t = elem (i, 0); \
+	      for (int j = 0; j < nc; j++) \
+		{ \
+		  retval.elem (i, j) = t; \
+		  if (j < nc - 1) \
+		    t OP elem (i, j+1); \
+		} \
+	    } \
+	} \
+      else \
+	{ \
+	  for (int j = 0; j < nc; j++) \
+	    { \
+	      ELT_TYPE t = elem (0, j); \
+	      for (int i = 0; i < nr; i++) \
+		{ \
+		  retval.elem (i, j) = t; \
+		  if (i < nr - 1) \
+		    t OP elem (i+1, j); \
+		} \
+	    } \
+	} \
+    } \
+ \
+  return retval
+
+#define MX_BASE_REDUCTION_OP(RET_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, \
+			     MT_RESULT) \
+ \
+  int nr = rows (); \
+  int nc = cols (); \
+ \
+  RET_TYPE retval; \
+ \
+  if (nr > 0 && nc > 0) \
+    { \
+      if ((nr == 1 && dim == -1) || dim == 1) \
+	{ \
+	  retval.resize (nr, 1); \
+	  for (int i = 0; i < nr; i++) \
+	    { \
+	      retval.elem (i, 0) = INIT_VAL; \
+	      for (int j = 0; j < nc; j++) \
+		{ \
+		  ROW_EXPR; \
+		} \
+	    } \
+	} \
+      else \
+	{ \
+	  retval.resize (1, nc); \
+	  for (int j = 0; j < nc; j++) \
+	    { \
+	      retval.elem (0, j) = INIT_VAL; \
+	      for (int i = 0; i < nr; i++) \
+		{ \
+		  COL_EXPR; \
+		} \
+	    } \
+	} \
+    } \
+  else \
+    { \
+      retval.resize (1, 1); \
+      retval.elem (0, 0) = MT_RESULT; \
+    } \
+ \
+  return retval
+
+#define MX_REDUCTION_OP_ROW_EXPR(OP) \
+  retval.elem (i, 0) OP elem (i, j)
+
+#define MX_REDUCTION_OP_COL_EXPR(OP) \
+  retval.elem (0, j) OP elem (i, j)
+
+#define MX_REDUCTION_OP(RET_TYPE, OP, INIT_VAL, MT_RESULT) \
+  MX_BASE_REDUCTION_OP (RET_TYPE, \
+			MX_REDUCTION_OP_ROW_EXPR (OP), \
+			MX_REDUCTION_OP_COL_EXPR (OP), \
+			INIT_VAL, MT_RESULT)
 #endif
 
 /*
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,7 @@
+2001-11-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* data.cc (DATA_REDUCTION): If no DIM arg, pass -1 to FCN.
+
 2001-11-09  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* oct-conf.h.in (FFTW_LIBS, LD_CXX, LIBOCT_PATHSEARCH,
--- a/src/data.cc
+++ b/src/data.cc
@@ -261,11 +261,11 @@
     { \
       octave_value arg = args(0); \
  \
-      int dim = (nargin == 1 ? 0 : args(1).int_value (true) - 1); \
+      int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
  \
       if (! error_state) \
 	{ \
-	  if (dim == 0 || dim == 1) \
+	  if (dim <= 1 && dim >= -1) \
 	    { \
 	      if (arg.is_real_type ()) \
 		{ \