Mercurial > hg > octave-avbm
changeset 8751:9f7ce4bf7650
optimize min/max functions
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 16 Feb 2009 08:52:00 +0100 |
parents | 8af4ba6b4216 |
children | 06b9903a029b |
files | liboctave/CNDArray.cc liboctave/ChangeLog liboctave/dNDArray.cc liboctave/fCNDArray.cc liboctave/fNDArray.cc liboctave/intNDArray.cc liboctave/lo-mappers.cc liboctave/mx-inlines.cc liboctave/oct-cmplx.h |
diffstat | 9 files changed, 257 insertions(+), 757 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CNDArray.cc +++ b/liboctave/CNDArray.cc @@ -692,191 +692,25 @@ ComplexNDArray ComplexNDArray::max (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return max (dummy_idx, dim); + return do_mx_minmax_op<ComplexNDArray> (*this, dim, mx_inline_max); } ComplexNDArray ComplexNDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return ComplexNDArray (); - - dr(dim) = 1; - - ComplexNDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - Complex tmp_max; - - double abs_max = octave_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_max = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_max)) - { - abs_max = std::abs(tmp_max); - break; - } - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - Complex tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - - double abs_tmp = std::abs (tmp); - - if (abs_tmp > abs_max) - { - idx_j = j; - tmp_max = tmp; - abs_max = abs_tmp; - } - } - - if (xisnan (tmp_max)) - { - result.elem (i) = Complex_NaN_result; - idx_arg.elem (i) = 0; - } - else - { - result.elem (i) = tmp_max; - idx_arg.elem (i) = idx_j; - } - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<ComplexNDArray> (*this, idx_arg, dim, mx_inline_max); } ComplexNDArray ComplexNDArray::min (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return min (dummy_idx, dim); + return do_mx_minmax_op<ComplexNDArray> (*this, dim, mx_inline_min); } ComplexNDArray ComplexNDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return ComplexNDArray (); - - dr(dim) = 1; - - ComplexNDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - Complex tmp_min; - - double abs_min = octave_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_min = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_min)) - { - abs_min = std::abs(tmp_min); - break; - } - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - Complex tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - - double abs_tmp = std::abs (tmp); - - if (abs_tmp < abs_min) - { - idx_j = j; - tmp_min = tmp; - abs_min = abs_tmp; - } - } - - if (xisnan (tmp_min)) - { - result.elem (i) = Complex_NaN_result; - idx_arg.elem (i) = 0; - } - else - { - result.elem (i) = tmp_min; - idx_arg.elem (i) = idx_j; - } - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<ComplexNDArray> (*this, idx_arg, dim, mx_inline_min); } NDArray
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,16 @@ +2009-02-16 Jaroslav Hajek <highegg@gmail.com> + + * oct-cmplx.h (operator <, operator >): New operators. + * mx-inlines.cc (OP_MINMAX_FCN, OP_MINMAX_FCN2, OP_MINMAX_FCNN): + New macros. + (mx_inline_min, mx_inline_max, do_mx_minmax_op): New overloaded + template functions. + * dNDArray (NDArray::min, NDArray::max): Use do_mx_minmax_op. + * fNDArray (FloatNDArray::min, FloatNDArray::max): Ditto. + * CNDArray (ComplexNDArray::min, ComplexNDArray::max): Ditto. + * fCNDArray (FloatComplexNDArray::min, FloatComplexNDArray::max): + Ditto. + 2009-02-16 Jaroslav Hajek <highegg@gmail.com> * chMatrix.cc (charMatrix::all, charMatrix::any): Use do_mx_red_op.
--- a/liboctave/dNDArray.cc +++ b/liboctave/dNDArray.cc @@ -731,157 +731,25 @@ NDArray NDArray::max (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return max (dummy_idx, dim); + return do_mx_minmax_op<NDArray> (*this, dim, mx_inline_max); } NDArray NDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return NDArray (); - - dr(dim) = 1; - - NDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - double tmp_max = octave_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_max = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_max)) - break; - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - double tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - else if (tmp > tmp_max) - { - idx_j = j; - tmp_max = tmp; - } - } - - result.elem (i) = tmp_max; - idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j; - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<NDArray> (*this, idx_arg, dim, mx_inline_max); } NDArray NDArray::min (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return min (dummy_idx, dim); + return do_mx_minmax_op<NDArray> (*this, dim, mx_inline_min); } NDArray NDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return NDArray (); - - dr(dim) = 1; - - NDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - double tmp_min = octave_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_min = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_min)) - break; - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - double tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - else if (tmp < tmp_min) - { - idx_j = j; - tmp_min = tmp; - } - } - - result.elem (i) = tmp_min; - idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j; - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<NDArray> (*this, idx_arg, dim, mx_inline_min); } NDArray
--- a/liboctave/fCNDArray.cc +++ b/liboctave/fCNDArray.cc @@ -687,191 +687,25 @@ FloatComplexNDArray FloatComplexNDArray::max (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return max (dummy_idx, dim); + return do_mx_minmax_op<FloatComplexNDArray> (*this, dim, mx_inline_max); } FloatComplexNDArray FloatComplexNDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return FloatComplexNDArray (); - - dr(dim) = 1; - - FloatComplexNDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - FloatComplex tmp_max; - - float abs_max = octave_Float_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_max = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_max)) - { - abs_max = std::abs(tmp_max); - break; - } - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - FloatComplex tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - - float abs_tmp = std::abs (tmp); - - if (abs_tmp > abs_max) - { - idx_j = j; - tmp_max = tmp; - abs_max = abs_tmp; - } - } - - if (xisnan (tmp_max)) - { - result.elem (i) = FloatComplex_NaN_result; - idx_arg.elem (i) = 0; - } - else - { - result.elem (i) = tmp_max; - idx_arg.elem (i) = idx_j; - } - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<FloatComplexNDArray> (*this, idx_arg, dim, mx_inline_max); } FloatComplexNDArray FloatComplexNDArray::min (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return min (dummy_idx, dim); + return do_mx_minmax_op<FloatComplexNDArray> (*this, dim, mx_inline_min); } FloatComplexNDArray FloatComplexNDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return FloatComplexNDArray (); - - dr(dim) = 1; - - FloatComplexNDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - FloatComplex tmp_min; - - float abs_min = octave_Float_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_min = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_min)) - { - abs_min = std::abs(tmp_min); - break; - } - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - FloatComplex tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - - float abs_tmp = std::abs (tmp); - - if (abs_tmp < abs_min) - { - idx_j = j; - tmp_min = tmp; - abs_min = abs_tmp; - } - } - - if (xisnan (tmp_min)) - { - result.elem (i) = FloatComplex_NaN_result; - idx_arg.elem (i) = 0; - } - else - { - result.elem (i) = tmp_min; - idx_arg.elem (i) = idx_j; - } - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<FloatComplexNDArray> (*this, idx_arg, dim, mx_inline_min); } FloatNDArray
--- a/liboctave/fNDArray.cc +++ b/liboctave/fNDArray.cc @@ -686,157 +686,25 @@ FloatNDArray FloatNDArray::max (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return max (dummy_idx, dim); + return do_mx_minmax_op<FloatNDArray> (*this, dim, mx_inline_max); } FloatNDArray FloatNDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return FloatNDArray (); - - dr(dim) = 1; - - FloatNDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - float tmp_max = octave_Float_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_max = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_max)) - break; - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - float tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - else if (tmp > tmp_max) - { - idx_j = j; - tmp_max = tmp; - } - } - - result.elem (i) = tmp_max; - idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j; - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<FloatNDArray> (*this, idx_arg, dim, mx_inline_max); } FloatNDArray FloatNDArray::min (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return min (dummy_idx, dim); + return do_mx_minmax_op<FloatNDArray> (*this, dim, mx_inline_min); } FloatNDArray FloatNDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = dims (); - dim_vector dr = dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return FloatNDArray (); - - dr(dim) = 1; - - FloatNDArray result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j; - - float tmp_min = octave_Float_NaN; - - for (idx_j = 0; idx_j < x_len; idx_j++) - { - tmp_min = elem (idx_j * x_stride + x_offset); - - if (! xisnan (tmp_min)) - break; - } - - for (octave_idx_type j = idx_j+1; j < x_len; j++) - { - float tmp = elem (j * x_stride + x_offset); - - if (xisnan (tmp)) - continue; - else if (tmp < tmp_min) - { - idx_j = j; - tmp_min = tmp; - } - } - - result.elem (i) = tmp_min; - idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j; - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<FloatNDArray> (*this, idx_arg, dim, mx_inline_min); } FloatNDArray
--- a/liboctave/intNDArray.cc +++ b/liboctave/intNDArray.cc @@ -213,140 +213,28 @@ intNDArray<T> intNDArray<T>::max (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return max (dummy_idx, dim); + return do_mx_minmax_op<intNDArray<T> > (*this, dim, mx_inline_max); } template <class T> intNDArray<T> intNDArray<T>::max (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = this->dims (); - dim_vector dr = this->dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return intNDArray<T> (); - - dr(dim) = 1; - - intNDArray<T> result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j = 0; - - T tmp_max = this->elem (x_offset); - - for (octave_idx_type j = 1; j < x_len; j++) - { - T tmp = this->elem (j * x_stride + x_offset); - - if (tmp > tmp_max) - { - idx_j = j; - tmp_max = tmp; - } - } - - result.elem (i) = tmp_max; - idx_arg.elem (i) = idx_j; - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<intNDArray<T> > (*this, idx_arg, dim, mx_inline_max); } template <class T> intNDArray<T> intNDArray<T>::min (int dim) const { - ArrayN<octave_idx_type> dummy_idx; - return min (dummy_idx, dim); + return do_mx_minmax_op<intNDArray<T> > (*this, dim, mx_inline_min); } template <class T> intNDArray<T> intNDArray<T>::min (ArrayN<octave_idx_type>& idx_arg, int dim) const { - dim_vector dv = this->dims (); - dim_vector dr = this->dims (); - - if (dv.numel () == 0 || dim > dv.length () || dim < 0) - return intNDArray<T> (); - - dr(dim) = 1; - - intNDArray<T> result (dr); - idx_arg.resize (dr); - - octave_idx_type x_stride = 1; - octave_idx_type x_len = dv(dim); - for (int i = 0; i < dim; i++) - x_stride *= dv(i); - - for (octave_idx_type i = 0; i < dr.numel (); i++) - { - octave_idx_type x_offset; - if (x_stride == 1) - x_offset = i * x_len; - else - { - octave_idx_type x_offset2 = 0; - x_offset = i; - while (x_offset >= x_stride) - { - x_offset -= x_stride; - x_offset2++; - } - x_offset += x_offset2 * x_stride * x_len; - } - - octave_idx_type idx_j = 0; - - T tmp_min = this->elem (x_offset); - - for (octave_idx_type j = 1; j < x_len; j++) - { - T tmp = this->elem (j * x_stride + x_offset); - - if (tmp < tmp_min) - { - idx_j = j; - tmp_min = tmp; - } - } - - result.elem (i) = tmp_min; - idx_arg.elem (i) = idx_j; - } - - result.chop_trailing_singletons (); - idx_arg.chop_trailing_singletons (); - - return result; + return do_mx_minmax_op<intNDArray<T> > (*this, idx_arg, dim, mx_inline_min); } /*
--- a/liboctave/lo-mappers.cc +++ b/liboctave/lo-mappers.cc @@ -222,44 +222,18 @@ // (double, double) -> double mappers. -// FIXME -- need to handle NA too? +// According to Matlab, is both args are NaN, the first one is returned. double xmin (double x, double y) { - if (x < y) - return x; - - if (y <= x) - return y; - - if (xisnan (x) && ! xisnan (y)) - return y; - else if (xisnan (y) && ! xisnan (x)) - return x; - else if (octave_is_NA (x) || octave_is_NA (y)) - return octave_NA; - else - return octave_NaN; + return xisnan (y) ? x : (x <= y ? x : y); } double xmax (double x, double y) { - if (x > y) - return x; - - if (y >= x) - return y; - - if (xisnan (x) && ! xisnan (y)) - return y; - else if (xisnan (y) && ! xisnan (x)) - return x; - else if (octave_is_NA (x) || octave_is_NA (y)) - return octave_NA; - else - return octave_NaN; + return xisnan (y) ? x : (x >= y ? x : y); } // complex -> complex mappers.
--- a/liboctave/mx-inlines.cc +++ b/liboctave/mx-inlines.cc @@ -467,6 +467,163 @@ OP_CUM_FCNN (mx_inline_cumsum) OP_CUM_FCNN (mx_inline_cumprod) +#define OP_MINMAX_FCN(F, OP) \ +template <class T> \ +void F (const T *v, T *r, octave_idx_type n) \ +{ \ + if (! n) return; \ + T tmp = v[0]; \ + octave_idx_type i = 1; \ + while (xisnan (tmp) && i < n) tmp = v[i++]; \ + for (i = 1; i < n; i++) \ + if (v[i] OP tmp) tmp = v[i]; \ + *r = tmp; \ +} \ +template <class T> \ +void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n) \ +{ \ + if (! n) return; \ + T tmp = v[0]; \ + octave_idx_type tmpi = 0; \ + octave_idx_type i = 1; \ + while (xisnan (tmp) && i < n) tmp = v[i++]; \ + for (i = 1; i < n; i++) \ + if (v[i] OP tmp) { tmp = v[i]; tmpi = i; }\ + *r = tmp; \ + *ri = tmpi; \ +} + +OP_MINMAX_FCN (mx_inline_min, <) +OP_MINMAX_FCN (mx_inline_max, >) + +// Row reductions will be slightly complicated. We will proceed with checks +// for NaNs until we detect that no row will yield a NaN, in which case we +// proceed to a faster code. + +#define OP_MINMAX_FCN2(F, OP) \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \ +{ \ + if (! n) return; \ + bool nan = false; \ + octave_idx_type j = 0; \ + for (octave_idx_type i = 0; i < m; i++) \ + { \ + r[i] = v[i]; \ + if (xisnan (v[i])) nan = true; \ + } \ + j++; v += m; \ + while (nan && j < n) \ + { \ + nan = false; \ + for (octave_idx_type i = 0; i < m; i++) \ + { \ + if (xisnan (v[i])) \ + nan = true; \ + else if (xisnan (r[i]) || v[i] OP r[i]) \ + r[i] = v[i]; \ + } \ + j++; v += m; \ + } \ + while (j < n) \ + { \ + for (octave_idx_type i = 0; i < m; i++) \ + if (v[i] OP r[i]) r[i] = v[i]; \ + j++; v += m; \ + } \ +} \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type *ri, \ + octave_idx_type m, octave_idx_type n) \ +{ \ + if (! n) return; \ + bool nan = false; \ + octave_idx_type j = 0; \ + for (octave_idx_type i = 0; i < m; i++) \ + { \ + r[i] = v[i]; ri[i] = j; \ + if (xisnan (v[i])) nan = true; \ + } \ + j++; v += m; \ + while (nan && j < n) \ + { \ + nan = false; \ + for (octave_idx_type i = 0; i < m; i++) \ + { \ + if (xisnan (v[i])) \ + nan = true; \ + else if (xisnan (r[i]) || v[i] OP r[i]) \ + { r[i] = v[i]; ri[i] = j; } \ + } \ + j++; v += m; \ + } \ + while (j < n) \ + { \ + for (octave_idx_type i = 0; i < m; i++) \ + if (v[i] OP r[i]) \ + { r[i] = v[i]; ri[i] = j; } \ + j++; v += m; \ + } \ +} + +OP_MINMAX_FCN2 (mx_inline_min, <) +OP_MINMAX_FCN2 (mx_inline_max, >) + +#define OP_MINMAX_FCNN(F) \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type l, \ + octave_idx_type n, octave_idx_type u) \ +{ \ + if (! n) return; \ + if (l == 1) \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, n); \ + v += n; r++; \ + } \ + } \ + else \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, l, n); \ + v += l*n; \ + r += l; \ + } \ + } \ +} \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type *ri, \ + octave_idx_type l, octave_idx_type n, octave_idx_type u) \ +{ \ + if (! n) return; \ + if (l == 1) \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, ri, n); \ + v += n; r++; ri++; \ + } \ + } \ + else \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, ri, l, n); \ + v += l*n; \ + r += l; ri += l; \ + } \ + } \ +} + +OP_MINMAX_FCNN (mx_inline_min) +OP_MINMAX_FCNN (mx_inline_max) + // Assistant function inline void @@ -542,6 +699,52 @@ return ret; } +template <class ArrayType> +inline ArrayType +do_mx_minmax_op (const ArrayType& src, int dim, + void (*mx_minmax_op) (const typename ArrayType::element_type *, + typename ArrayType::element_type *, + octave_idx_type, octave_idx_type, octave_idx_type)) +{ + octave_idx_type l, n, u; + dim_vector dims = src.dims (); + get_extent_triplet (dims, dim, l, n, u); + + // If the dimension is zero, we don't do anything. + if (dim < dims.length () && dims(dim) != 0) dims(dim) = 1; + dims.chop_trailing_singletons (); + + ArrayType ret (dims); + mx_minmax_op (src.data (), ret.fortran_vec (), l, n, u); + + return ret; +} + +template <class ArrayType> +inline ArrayType +do_mx_minmax_op (const ArrayType& src, Array<octave_idx_type>& idx, int dim, + void (*mx_minmax_op) (const typename ArrayType::element_type *, + typename ArrayType::element_type *, + octave_idx_type *, + octave_idx_type, octave_idx_type, octave_idx_type)) +{ + octave_idx_type l, n, u; + dim_vector dims = src.dims (); + get_extent_triplet (dims, dim, l, n, u); + + // If the dimension is zero, we don't do anything. + if (dim < dims.length () && dims(dim) != 0) dims(dim) = 1; + dims.chop_trailing_singletons (); + + ArrayType ret (dims); + if (idx.dims () != dims) idx = Array<octave_idx_type> (dims); + + mx_minmax_op (src.data (), ret.fortran_vec (), idx.fortran_vec (), + l, n, u); + + return ret; +} + // Avoid some code duplication. Maybe we should use templates. #define MX_CUMULATIVE_OP(RET_TYPE, ELT_TYPE, OP) \
--- a/liboctave/oct-cmplx.h +++ b/liboctave/oct-cmplx.h @@ -29,6 +29,24 @@ typedef std::complex<double> Complex; typedef std::complex<float> FloatComplex; +// The default comparison of complex number is to compare by abs, then by arg. +// FIXME: this could be speeded up significantly. +template <class T> +inline bool operator < (const std::complex<T>& a, + const std::complex<T>& b) +{ + T ax = std::abs (a), bx = std::abs (b); + return ax < bx || (ax == bx && std::arg (a) < std::arg (b)); +} + +template <class T> +inline bool operator > (const std::complex<T>& a, + const std::complex<T>& b) +{ + T ax = std::abs (a), bx = std::abs (b); + return ax > bx || (ax == bx && std::arg (a) > std::arg (b)); +} + #endif /*