Re: Implementation of symrcm

[David Bateman]

Michael Weitzel wrote:

> David Bateman schrieb:
> 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 ;-)

Hey, everyone here is pretty approachable, particularly if you start off
by sending code rather than complaints :-)

>> 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 ...

In the end I just tried to produce the same result as you previously
had, though if I could split the heap construction into two loops it
would probably be simpler. Though as I wasn't sure I understood all the
code I just ensured that the heap construction was the same.

In fact as I basically only used the form of the sparse matrix to
accelerate your code I imagine something similar might be done to your
original C code if you are still using it elsewhere...

> 
> 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);
> [...]

Its fine, though its more general to have a single return point from
Octave routines...

>>   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?!

Ok.

> 
>>       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?!
> 

empty_arg is defined in src/utils.cc and is just "(nr == 0 || nc == 0)".
Since nr is supposed to equal nc, I think with a little re-organization
we can get rid of empty_arg. Something like


  if (nr != nc)
    {
      gripe_square_matrix_required("symrcm");
      return retval;
    }

  if (nr == 0 && nc == 0)
    {
      retval = NDArray (dim_vector (1, 0));
      return retval;
    }


> @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.

Ok, I updated the copyright and committed the code..

D.