On Fri, Nov 9, 2012 at 2:45 PM, Frédéric Bastien
On Fri, Nov 9, 2012 at 3:32 AM, Nathaniel Smith
wrote: On Fri, Nov 9, 2012 at 6:18 AM, Nicolas SCHEFFER
wrote: Fred,
Thanks for the advice. The code will only affect the part in _dotblas.c where gemm is called. There's tons of check before that make sure both matrices are of ndim 2. We should check though if we can do these tricks in other parts of the function.
Otherwise: - I've built against ATLAS 3.10 - I'm happy to add a couple more test for C and F-contiguous. I'm not sure how to get the third type (strided), would you have an example?
def with_memory_order(a, order): assert order in ("C", "F", "discontig") assert a.ndim == 2 if order in ("C", "F"): return np.asarray(a, order=order) else: buf = np.empty((a.shape[0] * 2, a.shape[1] * 2), dtype=a.dtype) buf[::2, ::2] = a # This returns a view onto every other element of 'buf': result = buf[::2, ::2] assert not result.flags.c_contiguous and not result.flags.f_contiguous return result
The following test for instance checks integrity against multiarray.dot, which I believe is default when not compiled with BLAS. Dot is a hard function to test imho, so if anybody has ideas on what kind of test they'd like to see, please let me know.
If that's ok I might now be able to: - Check for more bugs, I need to dig a bit more in the gemm call, make sure everything is ok. - Create an issue on github and link to this discussion - Make a commit in a seperate branch - Move forward like that.
== import numpy as np from time import time from numpy.testing import assert_almost_equal
def test_dot_regression(): """ Test numpy dot by comparing with multiarray dot """ np.random.seed(7) a = np.random.randn(3, 3) b = np.random.randn(3, 2) c = np.random.randn(2, 3)
_dot = np.core.multiarray.dot
assert_almost_equal(np.dot(a, a), _dot(a, a)) assert_almost_equal(np.dot(b, c), _dot(b, c)) assert_almost_equal(np.dot(b.T, c.T), _dot(b.T, c.T))
assert_almost_equal(np.dot(a.T, a), _dot(a.T, a)) assert_almost_equal(np.dot(a, a.T), _dot(a, a.T)) assert_almost_equal(np.dot(a.T, a.T), _dot(a.T, a.T))
You should check that the result is C-contiguous in all cases too.
for a_order in ("C", "F", "discontig"): for b_order in ("C", "F", "discontig"): this_a = with_memory_order(a, a_order) this_b = with_memory_order(b, b_order) result = np.dot(this_a, this_b) assert_almost_equal(result, expected) assert result.flags.c_contiguous
You could also wrap the above in yet another loop to try a few different combinations of a and b matrices (perhaps after sticking the code into a utility function, like run_dot_tests(a, b, expected), so the indentation doesn't get out of hand ;-)). Then you can easily test some of the edge cases, like Nx1 matrices.
I agree that tests are needed the for Nx1 and variant cases. I saw blas error being raised with some blas version.
You also need to test with the output provided, so there is 3 loops:
for a_order in ("C", "F", "discontig", "neg"): for b_order in ("C", "F", "discontig", "neg"): for c_order in ("C", "F", "discontig", "neg"):
I also added the stride type "neg", I'm not sure if it is needed, but that is other corner cases. neg => result = buf[::-1, ::-1]
I just looked again at our code and there is another constrain: that the strides are multiple of the elemsize. Theano do not support not aligned array, but numpy does, so there is a need for test for this. You can make an unaligned array like this:
dtype = "b1,f4" a = numpy.empty(1e4, dtype=dtype)['f1']
I think you're mixing up two issues here... requiring that the strides and the itemsize match is part of the requirement for C- and F-contiguity, the example code you give here produces a non-contiguous array, so it's already checked for by the function we're talking about. However, there can be contiguous arrays that are not aligned, like: In [25]: a = np.empty(100, dtype="i1")[1:-3].view(np.float32) In [26]: a.flags.c_contiguous Out[26]: True In [27]: a.flags.aligned Out[27]: False I suspect the np.dot code that Nicolas is looking at already checks for the ALIGNED flag and makes a copy if necessary, but it would be good to have an array like this in the tests to be sure. -n
I just saw the strides problems that affect only some. I think that the best explaination is our code:
/* create appropriate strides for malformed matrices that are row or column * vectors, or empty matrices.
* In that case, the value of the stride does not really matter, but * some versions of BLAS insist that:
* - they are not smaller than the number of elements in the array, * - they are not 0. */
sx_0 = (Nx[0] > 1) ? Sx[0]/type_size : (Nx[1] + 1); sx_1 = (Nx[1] > 1) ? Sx[1]/type_size : (Nx[0] + 1);
sy_0 = (Ny[0] > 1) ? Sy[0]/type_size : (Ny[1] + 1); sy_1 = (Ny[1] > 1) ? Sy[1]/type_size : (Ny[0] + 1);
sz_0 = (Nz[0] > 1) ? Sz[0]/type_size : (Nz[1] + 1); sz_1 = (Nz[1] > 1) ? Sz[1]/type_size : (Nz[0] + 1);
So this ask for test with empty matrices too.
HTH
Fred
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion