Mercurial > hg > octave-max
changeset 9823:9b62f2d8de6d
improve r->c mapper strategy
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 17 Nov 2009 11:22:36 +0100 |
parents | 64270d3ad469 |
children | 6631c61a4a4e |
files | src/ChangeLog src/ov-flt-re-mat.cc src/ov-re-mat.cc |
diffstat | 3 files changed, 113 insertions(+), 20 deletions(-) [+] |
line wrap: on
line diff
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,10 @@ +2009-11-17 Jaroslav Hajek <highegg@gmail.com> + + * ov-re-mat.cc (do_rc_map): New static function. + (octave_matrix::map): Use it here. + * ov-flt-re-mat.cc (do_rc_map): New static function. + (octave_float_matrix::map): Use it here. + 2009-11-17 Jaroslav Hajek <highegg@gmail.com> * DLD-FUNCTIONS/balance.cc: Fix docs.
--- a/src/ov-flt-re-mat.cc +++ b/src/ov-flt-re-mat.cc @@ -680,6 +680,45 @@ return retval; } +// This uses a smarter strategy for doing the complex->real mappers. We +// allocate an array for a real result and keep filling it until a complex +// result is produced. +static octave_value +do_rc_map (const FloatNDArray& a, FloatComplex (&fcn) (float)) +{ + octave_idx_type n = a.numel (); + NoAlias<FloatNDArray> rr (a.dims ()); + + for (octave_idx_type i = 0; i < n; i++) + { + OCTAVE_QUIT; + + FloatComplex tmp = fcn (a(i)); + if (tmp.imag () == 0.0) + rr(i) = tmp.real (); + else + { + NoAlias<FloatComplexNDArray> rc (a.dims ()); + + for (octave_idx_type j = 0; j < i; j++) + rc(j) = rr(j); + + rc(i) = tmp; + + for (octave_idx_type j = i+1; j < n; j++) + { + OCTAVE_QUIT; + + rc(j) = fcn (a(i)); + } + + return new octave_float_complex_matrix (rc); + } + } + + return rr; +} + octave_value octave_float_matrix::map (unary_mapper_t umap) const { @@ -706,18 +745,22 @@ case umap_ ## UMAP: \ return octave_value (matrix.map<TYPE> (FCN)) - ARRAY_MAPPER (acos, FloatComplex, rc_acos); - ARRAY_MAPPER (acosh, FloatComplex, rc_acosh); +#define RC_ARRAY_MAPPER(UMAP, TYPE, FCN) \ + case umap_ ## UMAP: \ + return do_rc_map (matrix, FCN) + + RC_ARRAY_MAPPER (acos, FloatComplex, rc_acos); + RC_ARRAY_MAPPER (acosh, FloatComplex, rc_acosh); ARRAY_MAPPER (angle, float, ::arg); ARRAY_MAPPER (arg, float, ::arg); - ARRAY_MAPPER (asin, FloatComplex, rc_asin); + RC_ARRAY_MAPPER (asin, FloatComplex, rc_asin); ARRAY_MAPPER (asinh, float, ::asinhf); ARRAY_MAPPER (atan, float, ::atanf); - ARRAY_MAPPER (atanh, FloatComplex, rc_atanh); + RC_ARRAY_MAPPER (atanh, FloatComplex, rc_atanh); ARRAY_MAPPER (erf, float, ::erff); ARRAY_MAPPER (erfc, float, ::erfcf); ARRAY_MAPPER (gamma, float, xgamma); - ARRAY_MAPPER (lgamma, FloatComplex, rc_lgamma); + RC_ARRAY_MAPPER (lgamma, FloatComplex, rc_lgamma); ARRAY_MAPPER (ceil, float, ::ceilf); ARRAY_MAPPER (cos, float, ::cosf); ARRAY_MAPPER (cosh, float, ::coshf); @@ -725,16 +768,16 @@ ARRAY_MAPPER (expm1, float, ::expm1f); ARRAY_MAPPER (fix, float, ::fix); ARRAY_MAPPER (floor, float, ::floorf); - ARRAY_MAPPER (log, FloatComplex, rc_log); - ARRAY_MAPPER (log2, FloatComplex, rc_log2); - ARRAY_MAPPER (log10, FloatComplex, rc_log10); - ARRAY_MAPPER (log1p, FloatComplex, rc_log1p); + RC_ARRAY_MAPPER (log, FloatComplex, rc_log); + RC_ARRAY_MAPPER (log2, FloatComplex, rc_log2); + RC_ARRAY_MAPPER (log10, FloatComplex, rc_log10); + RC_ARRAY_MAPPER (log1p, FloatComplex, rc_log1p); ARRAY_MAPPER (round, float, xround); ARRAY_MAPPER (roundb, float, xroundb); ARRAY_MAPPER (signum, float, ::signum); ARRAY_MAPPER (sin, float, ::sinf); ARRAY_MAPPER (sinh, float, ::sinhf); - ARRAY_MAPPER (sqrt, FloatComplex, rc_sqrt); + RC_ARRAY_MAPPER (sqrt, FloatComplex, rc_sqrt); ARRAY_MAPPER (tan, float, ::tanf); ARRAY_MAPPER (tanh, float, ::tanhf); ARRAY_MAPPER (isna, bool, octave_is_NA);
--- a/src/ov-re-mat.cc +++ b/src/ov-re-mat.cc @@ -708,6 +708,45 @@ return retval; } +// This uses a smarter strategy for doing the complex->real mappers. We +// allocate an array for a real result and keep filling it until a complex +// result is produced. +static octave_value +do_rc_map (const NDArray& a, Complex (&fcn) (double)) +{ + octave_idx_type n = a.numel (); + NoAlias<NDArray> rr (a.dims ()); + + for (octave_idx_type i = 0; i < n; i++) + { + OCTAVE_QUIT; + + Complex tmp = fcn (a(i)); + if (tmp.imag () == 0.0) + rr(i) = tmp.real (); + else + { + NoAlias<ComplexNDArray> rc (a.dims ()); + + for (octave_idx_type j = 0; j < i; j++) + rc(j) = rr(j); + + rc(i) = tmp; + + for (octave_idx_type j = i+1; j < n; j++) + { + OCTAVE_QUIT; + + rc(j) = fcn (a(i)); + } + + return new octave_complex_matrix (rc); + } + } + + return rr; +} + octave_value octave_matrix::map (unary_mapper_t umap) const { @@ -734,18 +773,22 @@ case umap_ ## UMAP: \ return octave_value (matrix.map<TYPE> (FCN)) - ARRAY_MAPPER (acos, Complex, rc_acos); - ARRAY_MAPPER (acosh, Complex, rc_acosh); +#define RC_ARRAY_MAPPER(UMAP, TYPE, FCN) \ + case umap_ ## UMAP: \ + return do_rc_map (matrix, FCN) + + RC_ARRAY_MAPPER (acos, Complex, rc_acos); + RC_ARRAY_MAPPER (acosh, Complex, rc_acosh); ARRAY_MAPPER (angle, double, ::arg); ARRAY_MAPPER (arg, double, ::arg); - ARRAY_MAPPER (asin, Complex, rc_asin); + RC_ARRAY_MAPPER (asin, Complex, rc_asin); ARRAY_MAPPER (asinh, double, ::asinh); ARRAY_MAPPER (atan, double, ::atan); - ARRAY_MAPPER (atanh, Complex, rc_atanh); + RC_ARRAY_MAPPER (atanh, Complex, rc_atanh); ARRAY_MAPPER (erf, double, ::erf); ARRAY_MAPPER (erfc, double, ::erfc); ARRAY_MAPPER (gamma, double, xgamma); - ARRAY_MAPPER (lgamma, Complex, rc_lgamma); + RC_ARRAY_MAPPER (lgamma, Complex, rc_lgamma); ARRAY_MAPPER (ceil, double, ::ceil); ARRAY_MAPPER (cos, double, ::cos); ARRAY_MAPPER (cosh, double, ::cosh); @@ -753,16 +796,16 @@ ARRAY_MAPPER (expm1, double, ::expm1); ARRAY_MAPPER (fix, double, ::fix); ARRAY_MAPPER (floor, double, ::floor); - ARRAY_MAPPER (log, Complex, rc_log); - ARRAY_MAPPER (log2, Complex, rc_log2); - ARRAY_MAPPER (log10, Complex, rc_log10); - ARRAY_MAPPER (log1p, Complex, rc_log1p); + RC_ARRAY_MAPPER (log, Complex, rc_log); + RC_ARRAY_MAPPER (log2, Complex, rc_log2); + RC_ARRAY_MAPPER (log10, Complex, rc_log10); + RC_ARRAY_MAPPER (log1p, Complex, rc_log1p); ARRAY_MAPPER (round, double, xround); ARRAY_MAPPER (roundb, double, xroundb); ARRAY_MAPPER (signum, double, ::signum); ARRAY_MAPPER (sin, double, ::sin); ARRAY_MAPPER (sinh, double, ::sinh); - ARRAY_MAPPER (sqrt, Complex, rc_sqrt); + RC_ARRAY_MAPPER (sqrt, Complex, rc_sqrt); ARRAY_MAPPER (tan, double, ::tan); ARRAY_MAPPER (tanh, double, ::tanh); ARRAY_MAPPER (isna, bool, octave_is_NA);