Re: [PATCH 0 of 4] Implement basic sparse op diag and diag op sparse support.

[Jaroslav Hajek]

On Thu, Mar 12, 2009 at 3:39 PM, Jason Riedy <address@hidden> wrote:
> And Jaroslav Hajek writes:
>> I have applied all Jason's changes to my test repository, fixed the
>> diagonal * sparse problems, removed the checks for single-element
>> diagonal matrices and made a few other adjustments (such as not
>> introducing namespaces yet).
>
> Thanks!  I'll check on it when I have a sensible network connection
> tomorrow.  My work became hung up when I realized scalar double *
> complex diagonal was becoming double sparse * complex diagonal when
> performed explicitly, which is another reason why the single-element
> checks should be avoided.
>

Oh, this stopped working for you? It works without problems in my test repo.

> I became lost trying to add scalar double * complex diagonal[1].  ;)
>

That should be handled by op-cdm-s, I think. This only defines some
macros and includes op-dms-template.cc.
I don't understand where your problem was, but there doesn't seem to
be any. In any case, you shouldn't be able to create a 1x1 diagonal
matrix from the interpreter, it should be narrowed. You can force it
from a DLD function, but then you should know what you do, and there's
IMHO no need to make it still work as a scalar.

> At some point, the scalar->matrix constructors in liboctave should be
> marked as explicit.  That may flush out some surprising things that just
> happen to work...

Feel free to propose such changes. But maybe this should happen after
the anticipated API rewrite.

>> These changes do not yet modify the default conversion of diagonal
>> matrices - that would be the most incompatible change, probably, and
>> maybe should be further discussed.
>
> I still need to revisit your comments on the permutation matrices.
> Thanks for looking at all this!

I included those, too (that was more straightforward). The only thing
that is not handled is sparse + perm, but I don't really see a need
for that.

> With the basic sparse & diag (and sparse & permutation) ops, I don't see
> why either diag or sparse should fall back to sparse.  And the
> discussion of the "right" thing to do with diag .* sparse or sparse ./
> diag definitely is future work; there are at least two "right" things,
> depending on the use.
>
> Jason
>
> Footnotes:
> [1]  For one reason why that might be useful, consider
>        1 * diag ([1 + I() * Inf, 1])
> Converting the scalar 1 to a complex produces NaNs.  For further

This operation does not convert the multiplier 1 to complex. Try
omitting the "1 *" - it will do the same thing.
The problem is that i*Inf produces complex(NaN, Inf) because I() alone
gives complex(0,1).
It would probably be possible to optimize a+i*b to actually use
complex(a,b), but it's problematic because i may be overriden :(

> examples and reasons, see
>  W. Kahan and J.W. Thomas
>  Augmenting a Programming Language with Complex Arithmetic
>  http://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6127.html

That sounds interesting, thanks. We actually discussed having
"imaginary" values,  but it seemed too much work for too little
benefit.

> And yes, there's another whole discussion about Infs in complex
> arithmetic that leads back to a compromise made in the 70s. ;)
>

cheers


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