diagonal & permutation matrices

[Jaroslav Hajek]

hi all,

I'm working on patches allowing special treatment of diagonal and
permutation matrices in the Octave interpreter.

The motivation is that the natural syntax using diagonal & permutation
matrices, e.g.
D = norm (A, 'cols');
A /= D;
[Q,R,P] = qr (A); X = P' * R \ (Q' * B);
is currently formidable if performance is of any concern.
Multiplication & division by a permutation or diagonal matrix
is orders of magnitude slower that the "direct" counterpart
(subscripting by the permutation vector or using dmult).
And even though Matlab does the same, I think it's a bit shameful.

So my aim is to supply five octave_base_value subclasses:
octave_diag_matrix, octave_float_diag_matrix,
octave_complex_diag_matrix, octave_float_complex_diag_matrix,
octave_permutation_matrix
(it seems that a single permutation matrix subclass will be enough,
provided it will remember its numeric class)
and implement the necessary operations. Diagonal matrix support is
already present in liboctave, permutation matrix support needs to be
implemented (but I think it's no big deal, really).

I've run into two issues:

1. What to do when printing/saving? To honor Matlab compatibility, we
should display & save them as full matrices. However, it is certainly
much more informational to print them specially, and seems a little
bit weird for Octave to cripple its own data when saving variables.

2. If D is a diagonal matrix, the even D(i,i) = x with scalar i will
promote D to a full matrix. There does not seem to be an easy solution
avoiding this. (Just like even emtpy assignment to a range promotes it
to a matrix)

comments?

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