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

[Jaroslav Hajek]

On Tue, Mar 10, 2009 at 8:55 PM, Jason Riedy <address@hidden> wrote:
> And John W. Eaton writes:
>> At what point do we stop adding new features to 3.2?  How long should
>> we delay making the new release?
>
> It's up to you.  But these are not features that will be
> triggered in existing code, unless that code is relying on sparse
> * diag() or sparse * eye() becoming full.
>
> And Jaroslav Hajek writes:
>> Getting this into 3.2 would delay the release, which I think is not
>> good, especially when it seems that 3.0.x will have no more releases.
>
> Well, at least I'm flushing out bugs and adding missing tests.
>

That's nice. But I'm really glad you're interested in liboctave improvements.


>> 1. I'm not sure we need a separate "octave_impl" namespace. I was
>> hoping that just "octave" will suit well. Besides, the transition
>> should probably happen in an organized manner.
>
> I'll mangle the names in the next series rather than use a
> namespace.  I do not want these names escaping into the global
> symbol table without some protection.
>

Yes, that's what we do now. The result is, admittedly, a horrible mix-up
(think octave_idx_type vs. idx_vector, for example), but having a
namespace just in one header just doesn't seem right to me.


> BTW, you very much want a separate namespace to hold internal
> implementation details rather than splatting them all over a
> public namespace, and I'm assuming "octave" would be a great name
> for the public interface's namespace.

OK, maybe you're right. This can be discussed. What about having
octave::impl or octave::internal
instead? That seems better to me.

> And one strength of
> Octave's current design is that the changes can be incremental
> rather than made in one massive sweep.

I think there are strengths to both approaches. I'd certainly better
just distinguish between "old API" (what Octave has now, never much
official and never declared to be stable) and "new API" (in 3.4 or
4.0, whatever it will be, with namespaces and documented by doxygen or
something).

>> 2. I see you didn't reuse the existing mechanism Sparse-op-defs.h /
>> sparse-mx-ops / sparse-mk-ops.awk.
>
> As I mentioned in my commit message, the macros destroy compiler
> errors and debugging information.  Also, I cannot use
> sparse-mx-ops without hacking around element-wise multiplication
> and division.  I don't want to open that can right now.  The
> "right" result with structured matrices is open to many
> legitimate interpretations.
>
> So if I cannot use sparse-mx-ops without introducing changes that
> could destabilize other features, there's no point in sticking to
> the macrology.
>
> Also, touching Sparse-op-defs.h triggers a massive recompilation
> throughout liboctave *and* src.  Modifying Sparse-diag-op-defs.h
> only triggers rebuilding dSparse.o and CSparse.o along with
> relinking.  Very handy.
>

OK.

>> OK, the practice of putting executive code into macros is not really
>> nice, but I think you could have just create the necessary macros as
>> wrappers over the templates do_mul_dm_sm etc.
>
> That defeats my stated purpose.  If you insist, I'll do it, but
> it's a horrible idea.

No, it's fine, those are legitimate reasons. I just forgot that I kept
the support for diag & perm marices in mx-op-defs because a large part
was already there. I faced a similar issues and eventually I split up
mx-op-defs into two parts.

>> 3. You incorrectly implement the compound trans_mul and herm_mul
>> operations.
>
> Will fix, thanks!
>
>> 4. in the OPERATORS/ implementation, I don't think there's any need to
>> check for single-element diagonal and sparse matrices in operations.
>
> As I mentioned in my commit message, we need to ensure that the
> returned type is correct.  We cannot do so with the C++
> operator*.
>
>> A single element sparse [...] matrix will be narrowed to a scalar
>> automatically[...]
>
> No, not always.  Try it.
>> octave:1> sparse (1, 1, 1)
>> ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 1 [1e+02%])
>>
>>   (1, 1) ->  1
>> octave:2> sparse (1, 1, 1) * rand(2)
>> ans =
>>
>>    0.558860   0.409440
>>    0.942118   0.031385
>
> That's for Matlab(TM) compatibility.  A 1x1-dimension object does
> not *always* devolve into a scalar.
>
> And if it does, then I *still* need to fiddle with the types to
> ensure the correct type is returned.  Otherwise diag matrices
> don't act like full matrices.
>
>  sparse * scalar => sparse
>  sparse * full => full
>
>  sparse * 1x1 diag acting as a scalar => sparse
>  1x1 sparse acting as a scalar * diag acting as if full => DIAG
>
> We cannot make that distinction at the C++ operator* level.
>

OK, I see. Ugh. I don't really think there's any use for
sparse(scalar) * matrix, but it's OK to keep compatibility.
Still, at least the checking for single-element diag matrix can be
removed, because those will get properly narrowed, as well as dense
matrices (full-full or diag-full do not make such checks either).
That of course, rules out the reuse of op-dm-template. It was a nice
way to avoid code duplication, though.

>> Also, I don't understand
>> why you need to explicitly set matrix_type for the return value.
>
> I was trying to be clever and preserve the type as far as
> possible.  Will remove.

No, it's not a problem. I didn't realize that "mark_as_unsymmetric"
does something else than set the default. Probably my world is too
dense :)

cheers

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