sdot yields wrong results for odd sizes > 2^24 and all sizes > 2^29

[kohr-h]

When I compute a dot product of a vector of all ones with itself, the result should be the size of the vector. This starts to no longer hold at size 2^24 + 1 as the following program shows:

// file tmp.c
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "cblas.h"

int main(int argc, char **argv) {
    int exp = atoi(argv[1]);
    int offset = atoi(argv[2]);
    size_t size = (1 << exp) + offset;
    float *a = (float *) malloc(size * sizeof(float));
    assert(a != NULL);
    for (size_t i=0; i < size; ++i)
        a[i] = 1;

    float dot = cblas_sdot(size, a, 1, a, 1);
    printf("dot = %f\n", dot);
    printf("size = %zu = (1 << %i) + (%i)\n", size, exp, offset);

    return 0;
}

I get

$ ./tmp 24 0
dot = 16777216.000000
size = 16777216 = (1 << 24) + (0)

$ ./tmp 24 1
dot = 16777216.000000
size = 16777217 = (1 << 24) + (1)

$ ./tmp 24 2
dot = 16777218.000000
size = 16777218 = (1 << 24) + (2)

$ ./tmp 24 3
dot = 16777220.000000
size = 16777219 = (1 << 24) + (3)

$ ./tmp 24 4
dot = 16777220.000000
size = 16777220 = (1 << 24) + (4)

$ ./tmp 24 3
dot = 16777220.000000
size = 16777221 = (1 << 24) + (5)

and then, for sizes larger than 2^29 the result just stays constant at 2^29:

$ ./tmp 29 1
dot = 536870912.000000
size = 536870913 = (1 << 29) + (1)

$ ./tmp 29 2
dot = 536870912.000000
size = 536870914 = (1 << 29) + (2)

$ ./tmp 29 4
dot = 536870912.000000
size = 536870916 = (1 << 29) + (4)

$ ./tmp 30 -1
dot = 536870912.000000
size = 1073741823 = (1 << 30) + (-1)

$ ./tmp 30 0
dot = 536870912.000000
size = 1073741824 = (1 << 30) + (0)

I'm using the latest master develop version of OpenBLAS, and I statically link the library to not accidentally dynamically link an older version:

$ gcc -o tmp tmp.c [...]/OpenBLAS/build/lib/libopenblas.a -I[...]/OpenBLAS/build/include -std=c99
$ ldd tmp
        linux-vdso.so.1 (0x00007ffe2ed79000)
        libc.so.6 => /lib64/libc.so.6 (0x00007f58b2fdc000)
        /lib64/ld-linux-x86-64.so.2 (0x00007f58b338c000)

Some more info:

$ uname -a
Linux sabayon 4.11.0-sabayon #1 SMP Wed Jul 26 12:57:06 UTC 2017 x86_64 Intel(R) Core(TM) i7-4800MQ CPU @ 2.70GHz GenuineIntel GNU/Linux

See also this related Numpy issue

[brianborchers]

You're running into floating point round-off errors here. In single precision, IEEE floating point numbers have only 24 bits of precision, so the largest positive integer that can be represented precisely is 16777216=w^24.

You can verify that dot products of this size work fine in double precision.

[kohr-h]

Thanks for the explanation, that makes sense. Unfortunately that conclusion is a bit unsatisfactory since one would have to switch from float arrays to double arrays just to avoid this issue, giving away both speed and storage advantage of singles. This issue The second part of this issue would go away if the accumulator was double precision, at the expense of speed, of course.

Question: Is it completely out of scope of OpenBLAS to solve this issue by deciding upfront whether or not a double accumulator will be necessary or not?
I can imagine some simple heuristics, like computing the average of n << size samples and computing the expected dot value. But I know too little about floating-point arithmetic, and I also don't know the overall policy (try to return correct results vs. return predictably wrong results).
I do realize that the issue with the representation of the final result doesn't go away if the return value should be float, but the issue from 2^29 on seems to be more related to the summation itself, so that should be fixable.

[martin-frbg]

That suggestion does not sound too convincing to me - if you know exactly what loss of precision your particular use case can tolerate you may be better off coding your own case-specific trickery to save speed and memory ? There is an INTERFACE64 build option to switch all BLAS integers to long
but this will probably not save you if you want to preserve the speed and footprint of int32 (assuming it actually does what you need, which I have not checked now). I expect you will encounter this problem no matter which BLAS implementation you will try. Perhaps you can change your own code to use either single or double precision depending on size - that would keep all the advantages of floats where possible while still yielding a correct result for larger problems.

[adler-j]

A short note of interest from the related numpy issue: numpy/numpy#9852:

It seems this "bug" does not happen when using intelMKL, so as it stands OpenBLAS is in a way "inferior" to intelMKL, and downstream users get significantly different results depending on the backend used, perhaps this is something to keep in mind.

[martin-frbg]

What CPU are you seeing this with ?

[kohr-h]

That suggestion does not sound too convincing to me - if you know exactly what loss of precision your particular use case can tolerate you may be better off coding your own case-specific trickery to save speed and memory ?

No, that was just a lame suggestion, and I realize that I can do that easily myself.

There is an INTERFACE64 build option to switch all BLAS integers to long
but this will probably not save you if you want to preserve the speed and footprint of int32 (assuming it actually does what you need, which I have not checked now).

I don't see how this is related since we're not hitting any integer overflow issues yet. Of course, arrays larger than maxint are problematic as well, but that's a separate issue.

I expect you will encounter this problem no matter which BLAS implementation you will try.

When using Intel MKL with the same vanilla implementation as above, right, it produces the same results. However, in Numpy the dot product with MKL backend does not suffer from this issue, as opposed to the OpenBLAS backend. See numpy/numpy#9852.
I don't know what Numpy does on their end to fix this or if MKL exposes methods that let you specify precision of accumulator and result, for example. Maybe one the Numpy people can comment on that.
cuBLAS, for instance, has a dot function with extended capabilities for just that.

Perhaps you can change your own code to use either single or double precision depending on size - that would keep all the advantages of floats where possible while still yielding a correct result for larger problems.

Well, "you" is in this case Numpy, and I don't know about their plans. I as a user will have to resort to workarounds.

[kohr-h]

When using Intel MKL with the same vanilla implementation as above, right, it produces the same results.

Follow-up to this: the results are the same for the first half (sizes <= 2^28), but the issue with "constant result" from 2^28 + 1 on does not occur with MKL.

[martin-frbg]

I don't see how this is related since we're not hitting any integer overflow

That was just in the context of the odl issue you had linked (if I understood that one correctly)

As OpenBLAS uses different assembly kernels for ?dot on different CPUs, it would probably help to know which one(s) you are seeing this with.

[kohr-h]

That was just in the context of the odl issue you had linked (if I understood that one correctly)

Right, there are two entry points into this issue :-)

As OpenBLAS uses different assembly kernels for ?dot on different CPUs, it would probably help to know which one(s) you are seeing this with.

In the first post I specified everything that matters (I hope):

$ uname -a
Linux sabayon 4.11.0-sabayon #1 SMP Wed Jul 26 12:57:06 UTC 2017 x86_64 Intel(R) Core(TM) i7-4800MQ CPU @ 2.70GHz GenuineIntel GNU/Linux

and

    [...]
    float dot = cblas_sdot(size, a, 1, a, 1);
    [...]

[kohr-h]

Here's the last block of /proc/cpuinfo:

processor	: 7
vendor_id	: GenuineIntel
cpu family	: 6
model		: 94
model name	: Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz
stepping	: 3
microcode	: 0x39
cpu MHz		: 800.048
cache size	: 8192 KB
physical id	: 0
siblings	: 8
core id		: 3
cpu cores	: 4
apicid		: 7
initial apicid	: 7
fpu		: yes
fpu_exception	: yes
cpuid level	: 22
wp		: yes
flags		: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf tsc_known_freq pni pclmulqdq dtes64 monitor ds_cpl vmx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb intel_pt tpr_shadow vnmi flexpriority ept vpid fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt xsaveopt xsavec xgetbv1 xsaves dtherm ida arat pln pts hwp hwp_notify hwp_act_window hwp_epp
bugs		:
bogomips	: 8020.00
clflush size	: 64
cache_alignment	: 64
address sizes	: 39 bits physical, 48 bits virtual
power management:

[martin-frbg]

Thanks. I had indeed overlooked the uname info above, the cpuinfo is obviously from another system but both will be using the HASWELL target as far as OpenBLAS is concerned.

[kohr-h]

the cpuinfo is obviously from another system

Oops, that's me at my laptop yesterday and at my desktop machine today :-)

[brada4]

_dot is naive loop.
If you had even very random input you would get very different result if you sort input one or other way.
adding say in chunks could level out round-offs a bit (like hint cc heuristics), but i fail to see how that improves general cases.

Other options:
Make sdot with quad-precision sum register, all conversions are slow but accuracy impecable
convert en-masse to double and use ddot
neither is BLAS API anymore.

[kohr-h]

Note that there are actually two different issues here which I believe have different sources. The first is the lack of bits to represent the final result (the 2^24 - 2^28 case), which is baked into the API. The second one is the failure to accumulate beyond 2^28, which seems to be due to implementation rather than the API. So at least the second issue could be addressed without changing the API.

I still don't quite know what happens there. Can anyone solve the mystery for me?

[brada4]

Same limitation is present in reference implemetation.