compound operators - what next?

[Jaroslav Hajek]

Hello,

I've added optimized handling for expressions like A'*v and v*A' where
A is a sparse matrix and v a dense matrix or vector.
(useful, e.g., for some Krylov iterations).
I've also created a simple benchmark, that is attached.

On my laptop (Pentium M, gcc 4.2.1) I get with current Octave:

constructing sparse matrix
A'*v matrix: 1.294604 s
A'*v 100 vectors: 4.204328 s
(v'*A)' 100 vectors: 1.828778 s
v*A' matrix: 0.979149 s
v*A' 100 vectors: 4.621955 s

and with my repo:

constructing sparse matrix
A'*v matrix: 0.788072 s
A'*v 100 vectors: 0.841585 s
(v'*A)' 100 vectors: 1.388352 s
v*A' matrix: 0.305787 s
v*A' 100 vectors: 1.325876 s

on an AMD Opteron (2.6GHz, Intel C++), the results are:
constructing sparse matrix
A'*v matrix: 1.032092 s
A'*v 100 vectors: 3.318786 s
(v'*A)' 100 vectors: 1.376143 s
v*A' matrix: 0.632022 s
v*A' 100 vectors: 3.684349 s

and

A'*v matrix: 0.834610 s
A'*v 100 vectors: 0.830632 s
(v'*A)' 100 vectors: 1.122641 s
v*A' matrix: 0.603371 s
v*A' 100 vectors: 1.127208 s

thus especially the matrix-vector expressions are sped up dramatically
(which is obvious, since the transpose is avoided).
It is noteworthy that A'*v is faster than (v'*A)'.


So now, the following expressions are treated in an optimized way (A,B
dense matrices, S sparse):
A'*B A.'*B A*B' A*B.' (BLAS-GEMM)
A'*A A*A' (BLAS-SYRK, HERK)
S'*A A*S' (special code)
I'd like to ask anyone interested to try compiling, or review code,
see https://tw-math.de/highegg
What is needed to get these into mainstream? Shall I do a merge with main repo?

What next?
My top candidate is optimizing expressions like diag(v)*A because row
and column scaling tend to be quite common in my work (what about
yours?)
Since this will be more like a syntactic sugar (the above
optimizations really use different algorithms, or at least differently
optimized), I do not consider it that important; however, diag(v)*A is
IMHO more readable than dmult (v, A) (btw. what you do in Matlab?)
I see 3 ways open here:

1. implement compound operators diag_mul, mul_diag, diag_ldiv,
div_diag, so that the following expressions would get optimized:
diag(v)*A, A*diag(v), diag(v)\A, A/diag(v). Perhaps some checking will
be needed if `diag' refers to the intrinsic function.

2. this is suggested in Octave's old TODO: make v.*A work like
diag(v)*A. I don't like it much; but it's an option.

3. make `diag' and `eye' actually output special objects (i.e.
specialize octave_value for holding DiagMatrix). This would be the
most complex, since it would allow doing things like:
D = diag (v); % create scaling matrix
S1 = D*S;  % scale
S2 = S1 / (D+lambda*eye (n))
It would also partially eliminate the need for spdiag and speye,
because things like S+eye(n) or sparse(eye(n)) would work efficiently.
The question here is sensible diagonal matrix arithmetic vs. matlab
compatibility. For instance, should exp(D) for D diagonal only
exponentiate the diagonal? (this is sensible from a practical point of
view, but breaks Matlab compatibility).
Another question is whether this syntactic sugar is worth creating
another two octave_value classes (and thus increasing code bloat),
especially when most of the code will have to ensure compatibility.

I'll really appreciate anyone's comments on this issue.

Also, I'll welcome suggestions for more compound operators
optimizations, like, for instance mapping A+B*I directly to complex
(A,B)
(but is that any good?)

regards,

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

bench_sptransmul.m
Description: Text Data