complex comparison ops vs. Matlab

[Jaroslav Hajek]

hi all,

by this changeset:
http://hg.savannah.gnu.org/hgweb/octave/rev/7dafdb8b062f

I finished my project of refactoring dense binary operator
implementations in liboctave. The approach is similar to what I've
done before for reduction operations:
The innermost loops on pointer arrays are encapsulated as overloaded
template functions, providing good performance; appliers are provided
that take array arguments and return array, applying a given loop.
Using this machinery, the (MS|SM|MM|NDS|SND|NDND)_(BIN|CMP|BOOL)_OP(S)
macros in mx-op-defs were changed into trivial one-line wrappers, and
the *OPS1, *OPS2 alternatives are gone.

The total of mx-op-defs.h, mx-inlines.cc and MArray-defs.h is now 200
lines shorter than 3.2.x, and that also includes the new stuff for
in-place binary ops. It will be even more if also the VS_,SV_,and
VV_OPS are replaced (which was not immediately possible because
*Vector classes are not constructible from dim_vector).

The new implementation also probably creates less object code; as the
loop bodies are passed by pointer and thus can be shared amongst the
matrix and ndarray wrappers. I'm not sure whether the GNU linker
actually merges the instances, but it certainly has the possibility
(the last change cut down liboctave size by 330KB. It's still 400KB
larger than 3.2.x, but that's not fair comparison because there's a
lot of new stuff since 3.2.x).

The logical binary ops were sped up significantly by this change,
although the code is logically equal - probably GCC was not powerful
enough to optimize the previous code:

n = 4e3;
a = rand (n) < 0.5;
b = rand (n) < 0.5;
tic; a & b; toc
tic; a | b; toc
tic; a & !b; toc
tic; a | !b; toc
tic; !a | b; toc
tic; !a & b; toc

GCC 4.3.1, -O3 -march=native, Core 2 Duo @ 2.83 GHz:

3.2.x:

Elapsed time is 0.0800021 seconds.
Elapsed time is 0.0852289 seconds.
Elapsed time is 0.080343 seconds.
Elapsed time is 0.0791581 seconds.
Elapsed time is 0.0736232 seconds.
Elapsed time is 0.074126 seconds.

tip:

Elapsed time is 0.018801 seconds.
Elapsed time is 0.0197818 seconds.
Elapsed time is 0.026242 seconds.
Elapsed time is 0.020916 seconds.
Elapsed time is 0.0185721 seconds.
Elapsed time is 0.0178182 seconds.

One important point is that I changed the way comparing complex number works.
I think that Matlab compares only real parts of complex arrays when
involved in comparison operators; at least the previous Octave
implementation did so.
OTOH, Matlab defines a rigorous ordering on complex numbers as the
lexicographical ordering of [abs(z), arg(z)], which is used in sort,
min and max, and is shared with Octave.
I couldn't find the Matlab behavior documented anywhere, so I presume
it's a relic of the past and probably stemming from the fact that
Matlab is widely known to store the real and imaginary parts of arrays
separately - so it's the lazy do-nothing approach. However, it
simplifies the implementation significantly to consider one universal
ordering of complex numbers (defined in oct-cmplx.h) that gets
automatically used everywhere. As a result, Octave now uses the strict
abs/arg ordering of complex numbers for comparing complex arrays:

octave:1> 1+5i > 2+2i
ans =  1
octave:2> 1+i > 1-i
ans =  1

This means that various "obvious" invariants, such as

octave:3> a = sort (rand(100, 1) + i*rand(100,1)); all (a(1:99) <= a(2:100))
ans =  1
octave:4> a = rand(100, 1) + i*rand(100,1); all (a <= max (a(:)))
ans =  1

now hold for complex arrays as well (previously, this was not true).
Personally, I think this new behavior is more logical and consistent
and I vote for a change. I suppose few code depends on this and broken
scripts can be trivially fixed by wrapping operands in real () (which
I'd suggest anyway).

If I'm outvoted on this, I'd still like to keep that Matlab
compatibility mess out of liboctave, so I'll probably add a hack into
the interpreter to apply "real" to the operands first. It could also
be made conditional (for instance, keep the old behavior in the
braindamage mode), but I think John already said he doesn't like
global variables fundamentally influencing code flow.

Opinions? Comments?

regards

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz