Re: Implementation of symrcm

[Michael Weitzel]

David Bateman schrieb:
> John W. Eaton wrote:
>> On  7-May-2007, David Bateman wrote:
>>
>> | Michael Weitzel wrote:
>> | 
>> | > OCTAVE_LOCAL_BUFFER doesn't seem to work with type
>> | > bool - why?! I allocated a vector of bools manually.
>>
>> It fails because OCTAVE_LOCAL_BUFFER is currently just a macro that
>> uses std::vector<T> and the C++ standard allows std::vector<bool> to
>> be defined using a specialized bit vector.

I see. I wondered what's behind this macro. I imagined something like a
template class with an internal storage and a cast operator:

template< typename T > class localbuffer {
private:
        T * buf_;
public:
        inline localbuffer(unsigned int size) : buf_(new T[size]) {}
        inline ~localbuffer() { delete[] buf_; }
        inline operator T* () const { return buf_; }
};
#define OCTAVE_LOCAL_BUFFER(T,var,size) localbuffer<T>(var)(size)

>> | There are a few issues with this code.
>> | 
>> | [...]

I'm very happy with all your changes. My code was a bit clumsy since
this was also my first look on Octave's sources. I had a look at your
paper "Sparse matrix implementation in Octave" and didn't dare to ask
for help ;-)

>>
>> | John, if Michael is happy with the changed code, I'd suggest adding this
>> | to octave core as it is a core function of matlab.. Note that I've
>> | unconditionally included config.h so that the compile can get done
>> | outside of the source tree. This should be reincluded when the code is
>> | committed..
>>
>> I would also add:
>>
>>   In C++ there is no need to always write "struct CMK_node".  Instead,
>>   you can just use "CMK_Node" as the type once the struct is declared.

I'm aware ... this was a relict of the original C code.

>>
>>   Instead of passing pointers to scalar values, please use
>>   references, and instead of passing a struct by value, you probably
>>   want to pass it as a const reference if possible.  So for example,
>>   the function Q_enq should be
>>
>>     static inline void 
>>     Q_enq (CMK_Node *Q, octave_idx_type N, octave_idx_type& qh,
>>         octave_idx_type& qt, const CMK_Node& o)
>>     {        
>>       Q[qt] = o;
>>       qt = (qt+1) % N;
>>     }
>>
>>   and the call would be
>>
>>     Q_enq (Q, N, qh, qt, v);
> 
> I believe Michael's original code was in C and so he was probably trying
> to either do the minimum to the code or trying to keep it maintainable
> in the context of a dual use with another piece of code. If the first is
> the case, then the above should definitely be done. If the second is the
> case we should probably have a second file in C with the major part of
> this code in it.

you're right. This was originally C code. But I'm not a clear friend of
C++'s call-by-reference syntax, since in Q_enq(Q,N,&qh,&qh,v) it's
immediately clear that qh and qt might be modified. Q_enq(Q,N,qh,qh,v)
looks like call-by-value -- of course the actual implementation of Q_enq
is cleaner and easier to read if references are used and I'm absolutely
ok with this change.

> However, as I've made some pretty intrusive changes to his code, I
> suspect we should complete the conversion to C++ and have done so in the
> attached. I also attach a small benchmark script for Michael so that he
> can compare the difference to his original code. With this I get the
> speedups as shown in the table below relative to Michael's original
> code.

Your did absolutely convincing work and I'm sorry I didn't do a better
job. I had a hard time when searching for a C++ implementation to learn
from since most of my searches ended in wrappers for Fortran codes.
I haven't checked your changes to calc_degrees yet ... but if it still
sums the number of non-zero elements of row i of the symmetric matrix in
D[i] everything sould be okay. The results look fine ...

Maybe the initialization could be simplified (I guessed a lot here):
>   octave_value retval;
>   int nargin = args.length();
> 
>   if (nargin != 1)
>     {
>       print_usage ();
>       return retval;
>     }
> 
>   octave_value arg = args (0);
[...]
>   octave_idx_type nr = arg.rows ();
>   octave_idx_type nc = arg.columns ();
> 
>   if (nr == 0 && nc == 0)
>     {
>       retval = 1.;
this should be changed to

retval (0) = NDArray (dim_vector (1, 0));

Matlab returns a 0x0 matrix in this case. Since a 1xN matrix is returned
in case of a NxN matrix returning a 1x0 matrix here is maybe more
consistent?!

>       return retval;
>     }
> 
>   int arg_is_empty = empty_arg ("symrcm", nr, nc);
>   if (arg_is_empty < 0)
>     return retval;
>   // what is this good for?
>   if (arg_is_empty > 0)
>     return octave_value (Matrix (1, 1, 1.));

What does "arg_is_empty > 0" mean? I think it made no sense to return a
non-empty matrix here?!

@John:
> Also, once this is included with Octave, is there any objection to
> using the normal Octave copyright statement (still with Michael as the
> copyright holder, but GPL and stating
>
>  This file is part of Octave.
>
>  Octave is free software; you can redistribute [...]

Please change the license to GPL and add the statement. Since the rest
of Octave is also licensed under GPL it's okay for me.

regards
-- 
Michael Weitzel