Re: NDArray and Matrix prototypes

[Rik]

On 06/25/2012 08:32 AM, John W. Eaton wrote:
> On 23-Jun-2012, Rik wrote:
>
> | I've come across what looks like some unnecessary code in some of the
> | special functions.  The code for betainc is probably the easiest to show
> | off what is going on.  In lo-specfun.h there are multiple ways of calling
> | the function.
> | 
> | extern OCTAVE_API double betainc (double x, double a, double b);
> | extern OCTAVE_API Matrix betainc (double x, double a, const Matrix& b);
> | extern OCTAVE_API Matrix betainc (double x, const Matrix& a, double b);
> | extern OCTAVE_API Matrix betainc (double x, const Matrix& a, const Matrix& 
> b);
> | 
> | extern OCTAVE_API NDArray betainc (double x, double a, const NDArray& b);
> | extern OCTAVE_API NDArray betainc (double x, const NDArray& a, double b);
> | extern OCTAVE_API NDArray betainc (double x, const NDArray& a, const
> | NDArray& b);
> | 
> | extern OCTAVE_API Matrix betainc (const Matrix& x, double a, double b);
> | extern OCTAVE_API Matrix betainc (const Matrix& x, double a, const Matrix& 
> b);
> | extern OCTAVE_API Matrix betainc (const Matrix& x, const Matrix& a, double 
> b);
> | extern OCTAVE_API Matrix betainc (const Matrix& x, const Matrix& a, const
> | Matrix& b);
> | 
> | extern OCTAVE_API NDArray betainc (const NDArray& x, double a, double b);
> | extern OCTAVE_API NDArray betainc (const NDArray& x, double a, const
> | NDArray& b);
> | extern OCTAVE_API NDArray betainc (const NDArray& x, const NDArray& a,
> | double b);
> | extern OCTAVE_API NDArray betainc (const NDArray& x, const NDArray& a,
> | const NDArray& b);
> | 
> | In reality, there is only an actual algorithm for the function which takes
> | 3 double values and returns a double value.  The others merely loop across
> | rows and columns (for Matrix arguments) or use Fortran indexing (for
> | NDArrays) and call the betainc routine on each element of the array.
> | 
> | My question is whether it is really necessary to support both Matrix and
> | NDArray calling forms.
>
> That is probably left over from when there was a different class
> hierarchy for NDArray and Matrix objects.
>
> Maybe now these functions should all take Array<double> arguments and
> return Array<double> values?  And if needed, we could provide trivial
> wrappers that return Matrix and NDArray objects.
>
> Also, instead of writing something like
>
>   Array<double>
>   betainc (double x, double a, const Array<double>& b)
>   {
>     dim_vector dv = b.dims ();
>     octave_idx_type nel = dv.numel ();
>
>     Array<double> retval (dv);
>
>     for (octave_idx_type i = 0; i < nel; i++)
>       retval (i) = betainc (x, a, b(i));
>
>     return retval;
>   }
>
> we should probably use something like
>
>   Array<double>
>   betainc (double x, double a, const Array<double>& b)
>   {
>     dim_vector dv = b.dims ();
>     octave_idx_type nel = dv.numel ();
>
>     Array<double> retval (dv);
>
>     double *pretval = retval.fortran_vec ();
>
>     for (octave_idx_type i = 0; i < nel; i++)
>       *pretval++ = betainc (x, a, b(i));
>
>     return retval;
>   }
>
> instead to avoid checking the retval reference count at each
> iteration.  There's no need to change b(i) since b is const.
6/28/12

John,

I implemented the change for betainc and then converted betaincinv over to
the same format.  Benchmarking shows that avoiding the reference count
check saves 20 ms on 1e6 numbers.  Still, this was 1% and the difference
was statistically significant so I used the pointer to retval technique.

--Rik