Re: Octave and GSOC'2008

[Jaroslav Hajek]

On Thu, Feb 28, 2008 at 11:42 AM, John W. Eaton <address@hidden> wrote:
> On 28-Feb-2008, Jaroslav Hajek wrote:
>
>
> | That's exactly why I thought more about using just a special type for
>  | octave_value.
>  | For instance, a common construct is A'*B, where the value of A' is not 
> needed
>  | anymore after the statement. Given that evaluating A'*B is typically
>  | even slightly
>  | faster than A*B, transposing  and then multiplying is really a waste
>  | of time and memory.
>  |
>  | So I was thinking about
>  | DEFUN_DLD(times, ...
>  | )
>  |
>  | aA = args(0);
>  | if (aA.is_matrix_expr())
>  | {
>  |   eA = aA.matrix_expr_value();
>  |   switch (eA.expr_type())
>  |   {
>  |      case MatrixExpr::transpose:
>  |         /* specialized code for transpose, possibly DGEMM('T'...) */
>  |      case MatrixExpr::diag:
>  |         /* just scale rows */
>  |   }
>  | } else if (aA.is_matrix())
>  | {
>  |   A = aA.matrix_value()
>  | ....
>  |
>  |
>  | Once matrix_value() gets called, the lazy expression is evaluated and
>  | substituted by the real thing. The problem is that matrix_value() now
>  | affects the object for which it is called - I can imagine that causing
>  | a lot of problems elsewhere...
>
>  If you want to handle something like A'*B as a special case, then I
>  think you just want to be looking at the parse tree and recognizing
>  the expression tree that matches the pattern for a transposed
>  expression times another expression, then calling a special function
>  to evaluate that.  More generally, you might want to do some other
>  transformations, so having a general mechanism to do this rather than
>  doing it on an ad-hoc case-by-case basis would be better (and of
>  course, more work).
>
>  jwe
>
True, but I was thinking about just picking the most common cases,
(which A'*B certainly is), not providing lazy evaluation in general.
Moreover, It would be nice to allow to do this in user functions (at
least DEFUNs), which
do not normaly access the parse tree. When coding in m-files, I think
transpose is one of the
most frequent "temporary" operations used. I would really like Q'*A to
call DGEMM('T',...) instead of transposing and then call the slightly
slower DGEMM('N',..). (Actually, I don't think Octave calls BLAS for
matrix-matrix multiply at all, but you get the idea).
A similar thing is doing A \ B' etc. kron is also frequently used with
ones(..) to spread matrices, where multiplying is quite unnecessary
(admittedly, the speed impact is probably low in this case).

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