NDArray and Matrix prototypes

[John W. Eaton]

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.

jwe