diff options
Diffstat (limited to 'numpy')
| -rw-r--r-- | numpy/linalg/tests/test_linalg.py | 32 | ||||
| -rw-r--r-- | numpy/linalg/umath_linalg.c.src | 26 |
2 files changed, 52 insertions, 6 deletions
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py index aefb61b14..d40157904 100644 --- a/numpy/linalg/tests/test_linalg.py +++ b/numpy/linalg/tests/test_linalg.py @@ -105,6 +105,9 @@ SQUARE_CASES = [ LinalgCase("8x8", np.random.rand(8, 8), np.random.rand(8)), + LinalgCase("1x1", + np.random.rand(1, 1), + np.random.rand(1)), LinalgCase("nonarray", [[1, 2], [3, 4]], [2, 1]), @@ -150,6 +153,12 @@ NONSQUARE_CASES = [ LinalgCase("8x11", np.random.rand(8, 11), np.random.rand(11)), + LinalgCase("1x5", + np.random.rand(1, 5), + np.random.rand(5)), + LinalgCase("5x1", + np.random.rand(5, 1), + np.random.rand(1)), ] HERMITIAN_CASES = [ @@ -177,7 +186,10 @@ HERMITIAN_CASES = [ None), LinalgCase("hmatrix_a_and_b", matrix([[1., 2.], [2., 1.]]), - None) + None), + LinalgCase("hmatrix_1x1", + np.random.rand(1, 1), + None), ] @@ -247,6 +259,24 @@ def _stride_comb_iter(x): assert np.all(xi == x) yield xi, "stride_" + "_".join(["%+d" % j for j in repeats]) + # generate also zero strides if possible + if x.ndim >= 1 and x.shape[-1] == 1: + s = list(x.strides) + s[-1] = 0 + xi = np.lib.stride_tricks.as_strided(x, strides=s) + yield xi, "stride_xxx_0" + if x.ndim >= 2 and x.shape[-2] == 1: + s = list(x.strides) + s[-2] = 0 + xi = np.lib.stride_tricks.as_strided(x, strides=s) + yield xi, "stride_xxx_0_x" + if x.ndim >= 2 and x.shape[:-2] == (1, 1): + s = list(x.strides) + s[-1] = 0 + s[-2] = 0 + xi = np.lib.stride_tricks.as_strided(x, strides=s) + yield xi, "stride_xxx_0_0" + for src in (SQUARE_CASES, NONSQUARE_CASES, HERMITIAN_CASES, diff --git a/numpy/linalg/umath_linalg.c.src b/numpy/linalg/umath_linalg.c.src index 3c2b316e2..b38038fc2 100644 --- a/numpy/linalg/umath_linalg.c.src +++ b/numpy/linalg/umath_linalg.c.src @@ -812,24 +812,32 @@ linearize_@TYPE@_matrix(void *dst_in, @typ@ *dst = (@typ@ *) dst_in; if (dst) { - int i; + int i, j; @typ@* rv = dst; fortran_int columns = (fortran_int)data->columns; fortran_int column_strides = (fortran_int)(data->column_strides/sizeof(@typ@)); fortran_int one = 1; for (i=0; i< data->rows; i++) { - if (column_strides >= 0) { + if (column_strides > 0) { FNAME(@copy@)(&columns, (void*)src, &column_strides, (void*)dst, &one); } - else { + else if (column_strides < 0) { FNAME(@copy@)(&columns, (void*)((@typ@*)src + (columns-1)*column_strides), &column_strides, (void*)dst, &one); } + else { + /* Zero stride has undefined behavior in some BLAS + implementations (e.g. OSX Accelerate), so do it + manually */ + for (j = 0; j < columns; ++j) { + memcpy((@typ@*)dst + j, (@typ@*)src, sizeof(@typ@)); + } + } src += data->row_strides/sizeof(@typ@); dst += data->columns; } @@ -855,17 +863,25 @@ delinearize_@TYPE@_matrix(void *dst_in, (fortran_int)(data->column_strides/sizeof(@typ@)); fortran_int one = 1; for (i=0; i < data->rows; i++) { - if (column_strides >= 0) { + if (column_strides > 0) { FNAME(@copy@)(&columns, (void*)src, &one, (void*)dst, &column_strides); } - else { + else if (column_strides < 0) { FNAME(@copy@)(&columns, (void*)src, &one, (void*)((@typ@*)dst + (columns-1)*column_strides), &column_strides); } + else { + /* Zero stride has undefined behavior in some BLAS + implementations (e.g. OSX Accelerate), so do it + manually */ + if (columns > 0) { + memcpy((@typ@*)dst, (@typ@*)src + (columns-1), sizeof(@typ@)); + } + } src += data->columns; dst += data->row_strides/sizeof(@typ@); } |