ILP64 OpenBLAS gives different result from regular OpenBLAS

[Zha0q1]

ILP64 OpenBLAS built with

CXX="g++ -fPIC" CC="gcc -fPIC" make -j DYNAMIC_ARCH=1 DYNAMIC_OLDER=1 USE_OPENMP=1 INTERFACE64=1 BINARY=64 NO_SHARED=0 NO_LAPACK=0

#include <iostream>
#include <cblas.h>
int main ( int argc, char* argv[] ) {
    const long n = ((long)1 << 31) - 1;
    std::cout << n <<std::endl;
    float*  A = new float[n];
    float*  B = new float[n];
    float*  C = new float[1];
    for(long i =0; i <n; i++){
        A[i] = 1;
        B[i] = 1;

    }
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, 1, n,  1.0, A, n, B, 1, 0.0, C, 1);
    std::cout << *C <<std::endl;
    delete[] A;
    delete[] B;
    delete[] C;
    return 0;
}

Compiled with g++ gemm.cpp -I /usr/local/include/ -lopenblas

On ILP64 OpenBLAS machine:

2147483647
1.93274e+09

On vanilla OpenBLAS machine:

2147483647
2.14748e+09

I think I am linking correctly, because if I change n to (long)1 << 31) + 1, then
ILP64

2147483649
1.93274e+09

Vanilla

2147483649
 ** On entry to SGEMM  parameter number  5 had an illegal value
0

[brada4]

Is it same CPU core type?
vanilla result looks like pedantic addition, other might be FMA or other way around.
for completeness you could calculate same in the loop acc=0.0 ; for i in 1..n acc = acc+(a[i]*b[[i]) ;;

[martin-frbg]

Make sure that you use the right flavor of include directory and LD_LIBRARY_PATH for your tests as well. I can not reproduce your problem on Haswell-class hardware at least.

[Zha0q1]

Is it same CPU core type?
vanilla result looks like pedantic addition, other might be FMA or other way around.
for completeness you could calculate same in the loop acc=0.0 ; for i in 1..n acc = acc+(a[i]*b[[i]) ;;

Yes they are the same machine type.
I tried to add in a loop:

#include <iostream>
#include <cblas.h>
int main ( int argc, char* argv[] ) {
    const int64_t n = ((int64_t)1 << 31) - 1;
    std::cout << n <<std::endl;
    float*  A = new float[n];
    float*  B = new float[n];
    float*  C = new float[1];
    for(long i =0; i <n; i++){
        A[i] = 1;
        B[i] = 1;

    }
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, 1, n,  1.0, A, n, B, 1, 0.0, C, 1);
    std::cout << *C <<std::endl;

    float acc = 0.0;
    for (long i=0; i<n; i++)
        acc += A[i] * B[i];
    std::cout << acc << std::endl;

    delete[] A;
    delete[] B;
    delete[] C;
    return 0;
}

the results are:
ILP64 OpenBLAS

2147483647
1.93274e+09
1.67772e+07

Vanilla OpenBLAS

2147483647
2.14748e+09
1.67772e+07

Make sure that you use the right flavor of include directory and LD_LIBRARY_PATH for your tests as well. I can not reproduce your problem on Haswell-class hardware at least.

Do I need to add any flags? g++ gemm.cpp -I /usr/local/include/ -lopenblas. I checked at the executable is linking to the correct libopenblas.so

[brada4]

@Zha0q1 you can check library loaded by ld.so with ldd command
What the naive loop reaches is limit of "significand"
https://en.wikipedia.org/wiki/Floating-point_arithmetic#IEEE_754:_floating_point_in_modern_computers
i.e 2^24 (it rounds properly even near the end) i.e close to 1,6777216 ^ 10^7

setting acc to double will give result matching expected integer solution.
e.g lua programming language has no integers, it says double-s represent particular integer range exactly.

I am puzzled why IPL64 and "normal" versions do not use same blocking algorithm to get similar results. i.e if you add 64 elements in same loop then to accumulator you get more, like in "normal" example

[Zha0q1]

@Zha0q1 you can check library loaded by ld.so with ldd command
What the naive loop reaches is limit of "significand"
https://en.wikipedia.org/wiki/Floating-point_arithmetic#IEEE_754:_floating_point_in_modern_computers
i.e 2^24 (it rounds properly even near the end) i.e close to 1,6777216 ^ 10^7

setting acc to double will give result matching expected integer solution.
e.g lua programming language has no integers, it says double-s represent particular integer range exactly.

I am puzzled why IPL64 and "normal" versions do not use same blocking algorithm to get similar results. i.e if you add 64 elements in same loop then to accumulator you get more, like in "normal" example

After changing acc to double I got:

2147483647
1.93274e+09
2.14748e+09

Also ldd gives:

(nomkl) ubuntu@ip-172-31-10-124:~$ ldd /lib/x86_64-linux-gnu/libdl.so.2
	linux-vdso.so.1 (0x00007ffe471a3000)
	libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f7e96daa000)
	/lib64/ld-linux-x86-64.so.2 (0x00007f7e9739f000)

I am puzzled why IPL64 and "normal" versions do not use same blocking algorithm to get similar results. i.e if you add 64 elements in same loop then to accumulator you get more, like in "normal" example

Sorry I still don't quire understand. Is there a way to fix this?

[brada4]

Thats property of floating point numbers, depending on algorithm you may saturate significand or get the arihmetically correct result.
There is nothing to fix, any result between those is "correct". use type with enough significand if you intend to sum wide range of values, like double served in this case, or quad will serve to int64.

[Zha0q1]

Thats property of floating point numbers, depending on algorithm you may saturate significand or get the arihmetically correct result.
There is nothing to fix, any result between those is "correct". use type with enough significand if you intend to sum wide range of values, like double served in this case, or quad will serve to int64.

But when I use vanilla openblas it does give the arithmetically correct result on the same input. Is there a way to make ILP openblas consistent with that?

[brada4]

reference BLAS gives 16megs result, thats the absolute truth in this case
Using dgemm via scripting gets int32 biggest

R> a<-matrix(1,2**31-1,1)
R> crossprod(a)

(needs 32+GB RAM to get through)

[brada4]

The fix you submitted in numpy gets the aritmetic result, its c++ equivalent is using dgemm.

[Zha0q1]

The fix you submitted in numpy gets the aritmetic result, its c++ equivalent is using dgemm.

Yeah right, but with vanilla openblas even sgemm gives correct result

[brada4]

Yeah right, but with vanilla openblas even sgemm gives correct result

Yep, but I dont get clue how ;-)

[Zha0q1]

Make sure that you use the right flavor of include directory and LD_LIBRARY_PATH for your tests as well. I can not reproduce your problem on Haswell-class hardware at least.

Are you trying with ILP 64 build?

[brada4]

Actually sgemv clamps in same place but sdot gets a bit further thanks to blocking being implied.