diagonal matrices specializations

[Jaroslav Hajek]

hi all,

this is another proposed improvement in the interpreter - diagonal
matrix objects.

Rationale:
Diagonal matrices are common objects in linear algebra. They are
frequently used to perform scaling and unscaling, and are also
naturally yielded by the svd and eig factorizations. In Octave and
Matlab alike, however, direct use of these matrices suffers a
significant performance penalty, because they're treated as full.
So, while expressions like `diag(D) *A' are common in manuals and
examples, they're avoided in practice by using something ugly like `A
.* repmat (D, 1, columns(A))'. Octave has introduced `dmult' for this
purpose, however, this function is not much used because it is not
Matlab compatible and it is not much known even to Octave users. And
anyway, the whole point of having a language like Octave is to allow
user to type equations as directly as possible, with the interpreter
taking care of performance.

Summary of changes:
There is existing support for rectangular diagonal matrices in
liboctave. (Yeah, rectangular. Search me why, but after a while I
started to like the idea, it's not much additional work and it fits
more uses nicely, SVD for instance). Some modifications were made to
improve performance (e.g., inline element access) and extend
functionality.
In the interpreter code, 4 new octave_value classes are introduced (
real/complex diagonal matrix + float counterparts), with appropriate
operators registered and specialized in OPERATORS. Several builtins
were modified to take advantage of the new functionality (diag, eye,
inv) while others, due to their fine coding, benefit automatically
(svd, eig). As a result, some tests need to be modified.

Two issues are left open:

1. Any assignment into a diagonal matrix converts it to a full matrix,
even an assigment to a single diagonal element. This is an inherent
consequence of current Octave's type system, where you can specialize
assignment based on types, but at the time when indices are resolved,
the octave_base_value polymorphic instance is already fixed (see
ops.h, octave_base_value::numeric_assign etc). I don't think this is
much of an issue, as diagonal matrices are seldom manipulated by
elements, but still it's a defect.

2. There was a prior debate how the special diagonal matrices should
be displayed & saved. Saving to external formats is out of question,
of course, but displaying and Octave's native saving is debatable.
Displaying a diagonal matrix by printing out just the diagonal (after
a notification of diagonality) is no doubt more informative than
displaying it as full. Similarly, when Octave saves a matrix in its
native format, it can be useful to preserve the diagonality
information rather than forgetting it. OTOH, there's compatibility.
In the patch presented, neither option is exploited, i.e. the matrix
is always converted to full.
My personal preference is to alter the printing and leave the saving,
but I'd like to hear more opinions about this.

In very near future, I intend to contribute a similar improvement
concerning permutation matrices, another common type of object in
linear algebra.

The patch (>100K) is available for download here:

http://artax.karlin.mff.cuni.cz/~hajej2am/ulozna/octave/diag-mat.diff

cheers

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