[Numpy-discussion] additional dtype argument to numpy.dot() (Re: Numpy-discussion Digest, Vol 29, Issue 48)

[David Henderson]

Hi Paul, list:

Thanks for the reply.

numpy.sum()  does indeed have a dtype for the accumulator for the  
sum. numpy.sum() does not implement an inner (dot) product, just a  
straight summation.

The feature I'm requesting is to add a similar accumulator type  
argument for numpy.dot().

After tracing some of the code, the there is a generic dot product  
algorithm for the various types. If a BLAS library is available, the  
BLAS routines will supplant the library, but only for the BLAS  
routines that exist.

BLAS doesn't support long double, so only a subset of the different  
calls for different types call BLAS.

The place where this feature comes into play is iterative refinement  
of a linear algebra solution. To implement iterative refinement, an  
error term is calculated as a matrix dot product with a higher  
precision than the elements of the matrix.

The floating point hardware on the i386/PPC/Sparc is ready to go for  
this application, its just that BLAS doesn't implement it.

The numpy.dot() function is missing the feature to supply the type of  
the accumulator for the sum.

I'm looking for a way to prototype some routines to perform iterative  
refinement, and python/numpy looks like an ideal platform.

numpy has _almost_ the feature set needed to implement this. Adding  
the accumulator type argument to numpy.dot() is a needed feature.

code snippet exhibiting float128 functionality that already exists:
--------------------------
Python 2.6.1 (r261:67515, Feb 12 2009, 16:56:09)
[GCC 4.0.1 (Apple Computer, Inc. build 5367)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
 >>> import numpy
 >>> z=numpy.longdouble((1,2,3))
 >>> z
array([1.0, 2.0, 3.0], dtype=float128)
 >>> y=numpy.double((4,5,6))
 >>> y
array([ 4.,  5.,  6.])
 >>> numpy.dot(y,y)
77.0
 >>> numpy.dot(y,z)
32.0
 >>> numpy.dot(z,z)
14.0
---------------------------

I'm prepared to make the changes to numpy. It looks like they would  
entail modifying the source file "multiarrayobject.c". The existing  
code is:

/*NUMPY_API
   Numeric.innerproduct(a,v)
*/
static PyObject *
PyArray_InnerProduct(PyObject *op1, PyObject *op2)
{

If I make changes, can I get them inserted into the numpy source base  
(after appropriate testing??)

Thanks in advance,

David Henderson


On Feb 14, 2009, at 8:04 AM, numpy-discussion-request@scipy.org wrote:

> Fri, 13 Feb 2009 12:04:12 -0800, David Henderson wrote:
>> I'd like accumulate the summation in extended precision, "double" sum
>> for float  inputs, "long double" sum for "double" inputs.
>>
>> A way of doing this is to add an optional third argument to dot - to
>> specify the summation type.
>
> `dot` does matrix-matrix, matrix-vector, or vector-vector products.  
> These
> are usually implemented via calling the BLAS linear algebra  
> library, and
> AFAIK BLAS only has routines for type-homogenous arguments. So I doubt
> this can be implemented for `dot` in any way.
>
> On the other hand, `numpy.sum` already has a `dtype` argument that
> specifies the accumulator data type. Maybe you can use that?
>
> -- 
> Pauli Virtanen
> -----------------------------
>