# HG changeset patch # User jwe # Date 1005915382 0 # Node ID e78705239df5fec3a7e317fcd56b95b82145c2e9 # Parent 55c8eee153cb170b06792f9bf48c616a09622c3f [project @ 2001-11-16 12:56:20 by jwe] diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc --- 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 diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,19 @@ +2001-11-16 John W. Eaton + + * 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 * Range.cc (Range::nelem_internal): Special case ranges that must diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc --- 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 diff --git a/liboctave/mx-inlines.cc b/liboctave/mx-inlines.cc --- 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 /* diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,7 @@ +2001-11-16 John W. Eaton + + * data.cc (DATA_REDUCTION): If no DIM arg, pass -1 to FCN. + 2001-11-09 John W. Eaton * oct-conf.h.in (FFTW_LIBS, LD_CXX, LIBOCT_PATHSEARCH, diff --git a/src/data.cc b/src/data.cc --- 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 ()) \ { \