Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DGEMM regression on SkylakeX #1955

Closed
staticfloat opened this issue Jan 8, 2019 · 85 comments
Closed

DGEMM regression on SkylakeX #1955

staticfloat opened this issue Jan 8, 2019 · 85 comments
Milestone

Comments

@staticfloat
Copy link
Contributor

staticfloat commented Jan 8, 2019

It looks like 45fe8cb has created a regression in Julia's pinv() calculations on SkylakeX. In particular, creating a Hilbert matrix of size 1000 x 100 and asking for the pseudo-inverse now calculates the wrong thing:

using LinearAlgebra

function hilb(T::Type, m::Integer, n::Integer)
    a = Matrix{T}(undef, m, n)
    for i=1:n
        for j=1:m
            a[j,i]=one(T)/(i+j-one(T))
        end
    end
    return a
end
hilb(m::Integer, n::Integer) = hilb(Float64,m,n)

a = hilb(1000, 100)
apinv = pinv(a)

Including the SkylakeX kernel gives the following answer:

100×1000 Array{Float64,2}:
  2.57526e6   -2.33247e6    2.21848e6    2.19307e6   -4.13046e6   …  -4.71439e5   -6.80621e5   -6.56864e5   -8.6676e5    -3.86363e5
 -1.22338e11  -2.20992e11   2.36372e11  -1.14835e11  -9.1049e10       1.13475e10   8.51702e9    3.51379e9    1.54455e10   2.8167e8
  2.45922e11   3.06366e11  -3.45368e11   6.99788e10   1.34305e11     -2.0333e10   -1.72032e10  -6.35715e9   -3.15537e10   3.42079e9
 -1.98151e10  -5.04668e10   6.22131e9   -4.37235e10  -3.29137e9       2.4302e9     3.26304e9    2.4001e9     5.31362e9    1.42072e8
 -3.96966e11   1.59586e10  -1.05208e11  -3.27214e11   3.74498e11      5.9171e10    8.17795e10   7.30775e10   1.11105e11   3.58937e10
 -1.15417e11  -1.8089e10    3.02927e10  -7.71434e10   6.99771e10  …   1.55784e10   1.9823e10    1.69518e10   2.74274e10   8.50587e9
 -7.91383e11   3.19284e11  -5.46209e11  -6.27492e11   1.00493e12      1.26433e11   1.8369e11    1.69215e11   2.44706e11   8.48811e10
 -1.12133e11  -3.60367e10  -1.90203e8   -9.61073e10   7.40389e10      1.55288e10   2.07537e10   1.779e10     2.90644e10   7.91371e9
  4.19346e11  -2.37896e11   3.70455e11   3.44512e11  -6.0224e11      -6.9822e10   -1.03364e11  -9.66593e10  -1.36339e11  -4.94829e10
 -9.51913e9   -1.29776e11   1.82371e11   1.45922e10  -1.32419e11     -3.82507e9   -1.03485e10  -1.22032e10  -1.15361e10  -6.94914e9
  2.81647e12  -3.70604e11   9.33114e11   2.35652e12  -2.88098e12  …  -4.29963e11  -5.98574e11  -5.4097e11   -8.06547e11  -2.73127e11
  2.23042e11  -1.50478e11   1.87991e11   1.81127e11  -3.31279e11     -3.78415e10  -5.5537e10   -5.24234e10  -7.25081e10  -2.8155e10
  8.06723e11  -3.62461e11   6.16052e11   6.57412e11  -1.0719e12      -1.30794e11  -1.91569e11  -1.77425e11  -2.54474e11  -8.93792e10
  5.53103e10  -1.38396e10  -1.52335e10   3.90026e10  -4.70444e10     -8.49163e9   -1.06612e10  -9.72492e9   -1.3894e10   -6.12903e9
  9.43701e11  -5.29597e11   7.95135e11   7.22863e11  -1.31683e12     -1.54868e11  -2.28251e11  -2.12445e11  -3.0164e11   -1.08025e11
 -4.21561e11   1.74302e10  -1.05198e11  -3.6234e11    4.02288e11  …   6.38291e10   8.7825e10    7.88592e10   1.18726e11   3.97235e10
  ⋮                                                               ⋱   ⋮
  6.68906e10  -6.53051e9    3.21447e10   4.33472e10  -6.42927e10     -9.61282e9   -1.3713e10   -1.20454e10  -1.88794e10  -5.21562e9
 -4.3416e10    9.89228e9   -1.06994e10  -3.73798e10   4.63849e10  …   6.77796e9    9.31398e9    8.54109e9    1.23867e10   4.64291e9
 -3.36704e10  -2.19847e10   3.36634e10  -2.54336e10   4.56403e9       4.11065e9    4.50039e9    3.55012e9    6.46187e9    1.91442e9
  1.00752e11  -9.78859e10   1.5892e11    7.5809e10   -1.8612e11      -1.78414e10  -2.82659e10  -2.68784e10  -3.70109e10  -1.31011e10
 -6.24935e11   1.72017e11  -3.18826e11  -5.14607e11   7.21608e11      9.79871e10   1.39412e11   1.27482e11   1.86508e11   6.45552e10
 -3.56138e11   4.75922e10  -1.07088e11  -2.96353e11   3.60894e11      5.43011e10   7.52645e10   6.80322e10   1.01325e11   3.46671e10
  1.41477e11  -1.55157e11   2.42801e11   1.08982e11  -2.78805e11  …  -2.57232e10  -4.11459e10  -3.94458e10  -5.35713e10  -1.95093e10
 -1.71703e11   6.41864e9   -2.29926e10  -1.45472e11   1.56983e11      2.57486e10   3.4882e10    3.13021e10   4.71032e10   1.62219e10
 -1.57418e11  -2.50531e10   1.42196e10  -1.16675e11   1.07792e11      2.18507e10   2.86411e10   2.47383e10   3.9572e10    1.19948e10
 -5.09849e11   1.45783e11  -2.73511e11  -3.97073e11   5.85682e11      7.9147e10    1.13045e11   1.02973e11   1.5164e11    5.11554e10
 -2.01401e11   9.45476e10  -1.55585e11  -1.48556e11   2.6349e11       3.21347e10   4.71454e10   4.34308e10   6.28048e10   2.14634e10
  6.55703e11  -1.91764e11   3.65181e11   5.46763e11  -7.75768e11  …  -1.03493e11  -1.481e11    -1.35717e11  -1.97988e11  -6.84949e10
  4.39019e11  -7.20111e10   1.66746e11   3.76386e11  -4.6781e11      -6.78667e10  -9.50602e10  -8.63599e10  -1.27719e11  -4.38774e10
 -1.35314e11   1.43088e11  -2.17464e11  -8.70376e10   2.51535e11      2.36837e10   3.76675e10   3.57591e10   4.93057e10   1.72898e10
  2.95696e10   5.86712e10  -7.47343e10   2.73326e10   3.03537e10     -2.52465e9   -1.2501e9    -4.13464e7   -2.64329e9    3.99575e7
 -7.70318e10  -3.93945e10   4.41432e10  -8.97948e10   4.01424e10      1.12653e10   1.37315e10   1.20607e10   1.87392e10   7.00317e9

Excluding the SkylakeX kernel (e.g. reverting to 544b069) gives the result:

100×1000 Array{Float64,2}:
     112.527      -6192.3         1.06925e5   -8.28373e5    3.21394e6   -6.01292e6   …  -2.99287e5   -3.02032e5   -3.04795e5   -3.07576e5
   -6305.8            4.64899e5  -9.07773e6    7.54681e7   -3.07426e8    5.99356e8       3.28027e7    3.31027e7    3.34047e7    3.37085e7
       1.1309e5      -9.42656e6   1.9735e8    -1.71896e9    7.25526e9   -1.46068e10     -8.71604e8   -8.79551e8   -8.8755e8    -8.95596e8
      -9.32272e5      8.33527e7  -1.82785e9    1.64741e10  -7.1497e10    1.47819e11      9.57181e9    9.65882e9    9.74639e9    9.83447e9
       3.98657e6     -3.73868e8   8.48896e9   -7.86389e10   3.49436e11  -7.39605e11     -5.19571e10  -5.24279e10  -5.29016e10  -5.33781e10
      -8.8007e6       8.57783e8  -2.00715e10   1.90643e11  -8.66324e11   1.8764e12   …   1.44167e11   1.45468e11   1.46778e11   1.48094e11
       7.90418e6     -8.06081e8   1.95704e10  -1.91875e11   8.97794e11  -2.00535e12     -1.75621e11  -1.77197e11  -1.78783e11  -1.80377e11
       2.40961e6     -2.157e8     4.61896e9   -3.98835e10   1.62513e11  -3.0448e11      -1.76662e9   -1.79326e9   -1.82037e9   -1.84804e9
      -5.54279e6      5.75778e8  -1.41936e10   1.40992e11  -6.67618e11   1.50955e12      1.37796e11   1.39031e11   1.40274e11   1.41523e11
      -4.00904e6      3.9166e8   -9.13369e9    8.60798e10  -3.86441e11   8.21383e11      5.04267e10   5.08898e10   5.13561e10   5.18251e10
       1.63339e6     -1.8959e8    5.10176e9   -5.45604e10   2.76223e11  -6.69193e11  …  -8.18735e10  -8.25983e10  -8.33273e10  -8.40596e10
       4.57405e6     -4.75446e8   1.17311e10  -1.16639e11   5.52734e11  -1.25049e12     -1.12594e11  -1.13605e11  -1.14622e11  -1.15644e11
       3.29825e6     -3.27272e8   7.73331e9   -7.3721e10    3.34416e11  -7.18287e11     -4.43866e10  -4.47954e10  -4.52068e10  -4.56208e10
  -37289.1            2.26601e7  -9.77279e8    1.3636e10   -8.32615e10   2.3651e11       4.70926e10   4.75026e10   4.79148e10   4.83286e10
      -2.81123e6      3.03789e8  -7.75788e9    7.96192e10  -3.89208e11   9.11377e11      9.80715e10   9.89445e10   9.98226e10   1.00705e11
      -3.65969e6      3.80616e8  -9.39927e9    9.35385e10  -4.43642e11   1.00441e12  …   8.90641e10   8.98655e10   9.06717e10   9.14823e10
       ⋮                                                                 ⋮           ⋱
      -5.52909e5      6.12481e7  -1.61047e9    1.70918e10  -8.69055e10   2.14183e11      3.86899e10   3.90212e10   3.93541e10   3.96884e10
 -992608.0            1.08145e8  -2.79998e9    2.92745e10  -1.46579e11   3.54876e11  …   5.63476e10   5.68342e10   5.73232e10   5.78142e10
      -1.35946e6      1.47079e8  -3.78272e9    3.92898e10  -1.95373e11   4.69163e11      6.97768e10   7.03821e10   7.09905e10   7.16015e10
      -1.59876e6      1.72258e8  -4.41282e9    4.56535e10  -2.26067e11   5.40144e11      7.69398e10   7.76094e10   7.82824e10   7.89583e10
      -1.68423e6      1.80867e8  -4.61853e9    4.76281e10  -2.35036e11   5.59247e11      7.6759e10    7.74288e10   7.81023e10   7.87786e10
      -1.5861e6       1.69768e8  -4.32122e9    4.44178e10  -2.18431e11   5.17537e11      6.82601e10   6.88577e10   6.94585e10   7.0062e10
      -1.28093e6      1.36538e8  -3.46127e9    3.54304e10  -1.73446e11   4.08663e11  …   5.09162e10   5.13641e10   5.18145e10   5.2267e10
      -7.85551e5      8.29992e7  -2.08563e9    2.1155e10   -1.02525e11   2.38529e11      2.56914e10   2.59206e10   2.61511e10   2.63827e10
      -1.22319e5      1.1641e7   -2.59973e8    2.29043e9   -9.22435e9    1.59028e10     -5.87688e9   -5.92236e9   -5.96796e9   -6.01363e9
       6.34345e5     -6.95395e7   1.8113e9    -1.90533e10   9.60279e10  -2.3435e11      -4.02598e10  -4.06052e10  -4.09522e10  -4.13007e10
       1.37734e6     -1.49026e8   3.83373e9   -3.98355e10   1.98204e11  -4.76404e11     -7.24162e10  -7.30427e10  -7.36725e10  -7.43049e10
       1.94231e6     -2.09213e8   5.35875e9   -5.54398e10   2.74569e11  -6.56287e11  …  -9.49885e10  -9.58134e10  -9.66426e10  -9.74753e10
       2.11244e6     -2.26877e8   5.79482e9   -5.97807e10   2.95167e11  -7.02918e11     -9.83543e10  -9.92107e10  -1.00072e11  -1.00936e11
       1.59804e6     -1.71043e8   4.3541e9    -4.47652e10   2.20217e11  -5.22078e11     -6.99398e10  -7.0551e10   -7.11654e10  -7.17825e10
   23549.4           -1.60268e6   1.77704e7    5.98452e7   -1.59272e9    7.57578e9       6.25708e9    6.3078e9     6.35871e9    6.40975e9
      -3.07648e6      3.31117e8  -8.47524e9    8.76252e10  -4.33698e11   1.03595e12      1.49876e11   1.51177e11   1.52485e11   1.53798e11

Note that the pinv() definition is using SVD internally, so this is turning into an LAPACK.gesdd() call, which is itself giving very different answers, so this should be easy to reproduce locally by passing a Hilbert matrix of the above dimensions in through whichever interface you wish to dgesdd.

@ViralBShah
Copy link
Contributor

It would be good to have a point release with a fix for this soon, if possible.

@ViralBShah
Copy link
Contributor

We're happy to test a patch on our entire testsuite and report before the release is tagged, in case other issues are uncovered.

@brianborchers
Copy link

It's possible that there's nothing actually wrong with the pseudoinverse computed using the new version of OpenBLAS.

This particular A matrix is extremely badly conditioned, which means that even tiny perturbations in A can result in huge changes to the exact arithmetic pseudoinverse. In actual double precision floating point, small changes in the order of arithmetic operations could result in dramatic changes to the computed pseudoinverse. The pseudoinverse isn't actually "bad" unless it has the wrong spectrum or otherwise fails to satisfy the Moore-Penrose defining properties.

FWIW, I did the same calculation in MATLAB (using MKL as the BLAS) on my Skylake system and got

>> A=hilb(1000); A=A(:,1:100);
>> P=pinv(A);
>> P(1:4,1:4)

ans =

  Columns 1 through 3

     8.738690667311354e+01    -3.720494093014765e+03     4.893343437332960e+04
    -3.788867195428541e+03     2.164369520845519e+05    -3.224537869197214e+06
     5.178776101104618e+04    -3.350230685340255e+06     5.360812947292385e+07
    -3.151417707054240e+05     2.191509727500118e+07    -3.680906556637023e+08

  Column 4

    -2.790072926842434e+05
     1.977629060403547e+07
    -3.452351522072226e+08
     2.458863285828566e+09

This is quite different from either of the pseudoinverses reported from Julia using different versions of OpenBLAS.

Since the Skylake AVX512 kernels will perform the arithmetic in different order from the older AVX kernels, there is good reason to expect that computations won't be bit for bit identical. Since this problem is so badly conditioned, "not bit for bit" basically means "can be completely different."

@ViralBShah
Copy link
Contributor

While Hilbert matrices are certainly ill-conditioned, it appears some other tests in our Bunch Kaufman testsuite are failing too.
https://gist.github.com/vtjnash/c4f09f3019b335a690862134807da41c

Also, the testsuite using Hilbert matrices is intentionally using these:
https://github.com/JuliaLang/julia/blob/c9e5a6a4878d55b6b3b5503038f874ea1d77c26f/stdlib/LinearAlgebra/test/pinv.jl#L89

Digging deeper.

@brianborchers
Copy link

The tests coded into test_pinv look good to me and I would be concerned about a failure. The trick will be to find a more easily reproducible problem.

@ViralBShah
Copy link
Contributor

There are lots of other tests failing, it turns out. I am trying to find a simple testcase to reproduce.

@staticfloat
Copy link
Contributor Author

Yep; I just took the simplest piece of the test I could find that showed a large difference between the two. The actual test is ensuring the Moore-Penrose conditions are met, e.g. norm(a*apinv*a-a)/norm(a) ≈ 0 (with respect to a practical numerical tolerance), and in this case the result we're getting is roughly 15 orders of magnitude off; the result is 8.5e11, when it is typically 2.9e-5. I would be surprised if your MATLAB pinv() gave a poor result for this norm test.

@martin-frbg
Copy link
Collaborator

Original issue was #1643 , too bad this did not get caught in 0.3.4 but I still do not have the hardware.

@brada4
Copy link
Contributor

brada4 commented Jan 9, 2019

@staticfloat is it possible to drill problem down to copy vs kernel routines?

@martin-frbg
Copy link
Collaborator

Yes, did you actually bisect to 45fe8cb or did you "only" jump to the last version that preceded any SkylakeX DGEMM changes ? PR #1793 was merged a few days later, and I notice it brought a few
"small changes" to the initial dgemm kernel as well as the addition of n/tcopy and dgemm_beta kernels.

@ViralBShah
Copy link
Contributor

@Keno can we give @martin-frbg access to our Skylake machine - Antarctic?

@staticfloat
Copy link
Contributor Author

I did actually bisect down to 45fe8cb; every version git bisect tested after that version failed. (I did not test all of them of course, but I tested 6-7 of them as git bisect worked its way backward from v0.3.5 to v0.3.3).

We can give you a login to a SkylakeX machine so you can investigate further. It is possible it is the copy routines and not the DGEMM routine, that is true. I thought that unlikely, but it is possible. My patch that fixed this issue for us in Julia comments out both the DGEMM kernel and the copy kernels; I'll try to see if I can narrow it down to just one kernel that is the issue.

@martin-frbg
Copy link
Collaborator

Thanks. I am not sure how successful I would be at debugging AVX512 code, but one simple thing to try would seem to be to use the gemm n/t copy_4 C codes instead of their _8 counterparts for DGEMMINCOPY/DGEMMITCOPY - this is how it is done on Haswell (which uses a "4x8" dgemm kernel as well) so there is at least a small chance that this could be an oversight from fenrus75's experiments with different kernel layouts.

@staticfloat
Copy link
Contributor Author

staticfloat commented Jan 9, 2019

Okay, it looks like if I disable these three lines by adding leading #'s, things work properly:

#DGEMMKERNEL    =  dgemm_kernel_4x8_skylakex.c

#DGEMMINCOPY    =  dgemm_ncopy_8_skylakex.c
#DGEMMITCOPY    =  dgemm_tcopy_8_skylakex.c

@brada4
Copy link
Contributor

brada4 commented Jan 9, 2019

You could copy respective options from known to work KERNEL.HASWEL

@martin-frbg
Copy link
Collaborator

@staticfloat if all three need to be commented this is pretty much equivalent to reverting 45fe8cb ... so does changing only the DGEMMINCOPY/DGEMMITCOPY to their _4 versions not help ?

@martin-frbg martin-frbg added this to the 0.3.6 milestone Jan 9, 2019
@staticfloat
Copy link
Contributor Author

staticfloat commented Jan 10, 2019

Setting this works okay:

#DGEMMKERNEL    =  dgemm_kernel_4x8_skylakex.c

DGEMMINCOPY    =  ../generic/gemm_ncopy_4.c
DGEMMITCOPY    =  ../generic/gemm_tcopy_4.c
DGEMMONCOPY    =  dgemm_ncopy_8_skylakex.c
DGEMMOTCOPY    =  dgemm_tcopy_8_skylakex.c

@martin-frbg
Copy link
Collaborator

So DGEMMKERNEL is definitely broken and needs to be commented out ? Will probably need help from @fenrus75 who wrote this kernel.

@fenrus75
Copy link
Contributor

back to test case.. is there anything I can easily test this with?

@fenrus75
Copy link
Contributor

(also @martin-frbg send me an email and we'll talk about skx hardware)

@fenrus75
Copy link
Contributor

ping? is there any way to reproduce a failure ?

@staticfloat
Copy link
Contributor Author

@fenrus75 do you need access to a SkylakeX machine? I can get you SSH access on a Linux SkylakeX machine if that would help.

@fenrus75
Copy link
Contributor

fenrus75 commented Jan 24, 2019 via email

@staticfloat
Copy link
Contributor Author

I work for Intel... in my cube alone I have 3 skylakeX's ;-)

Hah! I should have assumed from the Portland, OR in your profile. ;)

I'm going to need to call in the big linear algebra guns on this one; @andreasnoack could I ask you to help me distill this down into a simpler form (e.g. a single ccall to a BLAS or LAPACK function) for @fenrus75? The Julia code above should distill down into a single LAPACK function that then calls BLAS, but I would prefer to not get something wrong by making some assumption about the LinAlg internals.

@martin-frbg
Copy link
Collaborator

Yes, but the Makefile (and CMakeLists.txt) needs to make sure that your new test only gets included when LAPACK/LAPACKE is available (there is a precedent for this - test_potrs).

TiborGY added a commit to TiborGY/OpenBLAS that referenced this issue Jul 7, 2019
* With the Intel compiler on Linux, prefer ifort for the final link step 

icc has known problems with mixed-language builds that ifort can handle just fine. Fixes OpenMathLib#1956

* Rename operands to put lda on the input/output constraint list

* Fix wrong constraints in inline assembly

for OpenMathLib#2009

* Fix inline assembly constraints

rework indices to allow marking argument lda4 as input and output. For OpenMathLib#2009

* Fix inline assembly constraints

rework indices to allow marking argument lda as input and output.

* Fix inline assembly constraints

* Fix inline assembly constraints

* Fix inline assembly constraints in Bulldozer TRSM kernels

rework indices to allow marking i,as and bs as both input and output (marked operand n1 as well for simplicity). For OpenMathLib#2009

* Correct range_n limiting

same bug as seen in OpenMathLib#1388, somehow missed in corresponding PR OpenMathLib#1389

* Allow multithreading TRMV again

revert workaround introduced for issue OpenMathLib#1332 as the actual cause appears to be my incorrect fix from OpenMathLib#1262 (see OpenMathLib#1388)

* Fix error introduced during cleanup

* Reduce list of kernels in the dynamic arch build

to make compilation complete reliably within the 1h limit again

* init

* move fix to right place

* Fix missing -c option in AVX512 test

* Fix AVX512 test always returning false due to missing compiler option

* Make x86_32 imply NO_AVX2, NO_AVX512 in addition to NO_AVX

fixes OpenMathLib#2033

* Keep xcode8.3 for osx BINARY=32 build

as xcode10 deprecated i386

* Make sure that AVX512 is disabled in 32bit builds

for OpenMathLib#2033

* Improve handling of NO_STATIC and NO_SHARED

to avoid surprises from defining either as zero. Fixes OpenMathLib#2035 by addressing some concerns from OpenMathLib#1422

* init

* address warning introed with OpenMathLib#1814 et al

* Restore locking optimizations for OpenMP case

restore another accidentally dropped part of OpenMathLib#1468 that was missed in OpenMathLib#2004 to address performance regression reported in OpenMathLib#1461

* HiSilicon tsv110 CPUs optimization branch

add HiSilicon tsv110 CPUs  optimization branch

* add TARGET support for  HiSilicon tsv110 CPUs

* add TARGET support for HiSilicon tsv110 CPUs

* add TARGET support for HiSilicon tsv110 CPUs

* Fix module definition conflicts between LAPACK and ReLAPACK

for OpenMathLib#2043

* Do not compile in AVX512 check if AVX support is disabled

xgetbv is function depends on NO_AVX being undefined - we could change that too, but that combo is unlikely to work anyway

* ctest.c : add __POWERPC__ for PowerMac

* Fix crash in sgemm SSE/nano kernel on x86_64

Fix bug OpenMathLib#2047.

Signed-off-by: Celelibi <celelibi@gmail.com>

* param.h : enable defines for PPC970 on DarwinOS

fixes:
gemm.c: In function 'sgemm_':
../common_param.h:981:18: error: 'SGEMM_DEFAULT_P' undeclared (first use in this function)
 #define SGEMM_P  SGEMM_DEFAULT_P
                  ^

* common_power.h: force DCBT_ARG 0 on PPC970 Darwin

without this, we see
../kernel/power/gemv_n.S:427:Parameter syntax error
and many more similar entries

that relates to this assembly command
dcbt 8, r24, r18

this change makes the DCBT_ARG = 0
and openblas builds through to completion on PowerMac 970
Tests pass

* Make TARGET=GENERIC compatible with DYNAMIC_ARCH=1

for issue OpenMathLib#2048

* make DYNAMIC_ARCH=1 package work on TSV110.

* make DYNAMIC_ARCH=1 package work on TSV110

* Add Intel Denverton

for OpenMathLib#2048

* Add Intel Denverton

* Change 64-bit detection as explained in OpenMathLib#2056

* Trivial typo fix

as suggested in OpenMathLib#2022

* Disable the AVX512 DGEMM kernel (again)

Due to as yet unresolved errors seen in OpenMathLib#1955 and OpenMathLib#2029

* Use POSIX getenv on Cygwin

The Windows-native GetEnvironmentVariable cannot be relied on, as
Cygwin does not always copy environment variables set through Cygwin
to the Windows environment block, particularly after fork().

* Fix for OpenMathLib#2063: The DllMain used in Cygwin did not run the thread memory
pool cleanup upon THREAD_DETACH which is needed when compiled with
USE_TLS=1.

* Also call CloseHandle on each thread, as well as on the event so as to not leak thread handles.

* AIX asm syntax changes needed for shared object creation

* power9 makefile. dgemm based on power8 kernel with following changes : 32x unrolled 16x4 kernel and 8x4 kernel using (lxv stxv butterfly rank1 update). improvement from 17 to 22-23gflops. dtrmm cases were added into dgemm itself

* Expose CBLAS interfaces for I?MIN and I?MAX

* Build CBLAS interfaces for I?MIN and I?MAX

* Add declarations for ?sum and cblas_?sum

* Add interface for ?sum (derived from ?asum)

* Add ?sum

* Add implementations of ssum/dsum and csum/zsum

as trivial copies of asum/zsasum with the fabs calls replaced by fmov to preserve code structure

* Add ARM implementations of ?sum

(trivial copies of the respective ?asum with the fabs calls removed)

* Add ARM64 implementations of ?sum

as trivial copies of the respective ?asum kernels with the fabs calls removed

* Add ia64 implementation of ?sum

as trivial copy of asum with the fabs calls removed

* Add MIPS implementation of ?sum

as trivial copy of ?asum with the fabs calls removed

* Add MIPS64 implementation of ?sum

as trivial copy of ?asum with the fabs replaced by mov to preserve code structure

* Add POWER implementation of ?sum

as trivial copy of ?asum with the fabs replaced by fmr to preserve code structure

* Add SPARC implementation of ?sum

as trivial copy of ?asum with the fabs replaced by fmov to preserve code structure

* Add x86 implementation of ?sum

as trivial copy of ?asum with the fabs calls removed

* Add x86_64 implementation of ?sum

as trivial copy of ?asum with the fabs calls removed

* Add ZARCH implementation of ?sum

as trivial copies of the respective ?asum kernels with the ABS and vflpsb calls removed

* Detect 32bit environment on 64bit ARM hardware

for OpenMathLib#2056, using same approach as OpenMathLib#2058

* Add cmake defaults for ?sum kernels

* Add ?sum

* Add ?sum definitions for generic kernel

* Add declarations for ?sum

* Add -lm and disable EXPRECISION support on *BSD

fixes OpenMathLib#2075

* Add in runtime CPU detection for POWER.

* snprintf define consolidated to common.h

* Support INTERFACE64=1

* Add support for INTERFACE64 and fix XERBLA calls

1. Replaced all instances of "int" with "blasint"
2. Added string length as "hidden" third parameter in calls to fortran XERBLA

* Correct length of name string in xerbla call

* Avoid out-of-bounds accesses in LAPACK EIG tests

see Reference-LAPACK/lapack#333

* Correct INFO=4 condition

* Disable reallocation of work array in xSYTRF

as it appears to cause memory management problems (seen in the LAPACK tests)

* Disable repeated recursion on Ab_BR in ReLAPACK xGBTRF

due to crashes in LAPACK tests

* sgemm/strmm

* Update Changelog with changes from 0.3.6

* Increment version to 0.3.7.dev

* Increment version to 0.3.7.dev

* Misc. typo fixes

Found via `codespell -q 3 -w -L ith,als,dum,nd,amin,nto,wis,ba -S ./relapack,./kernel,./lapack-netlib`

* Correct argument of CPU_ISSET for glibc <2.5

fixes OpenMathLib#2104

* conflict resolve

* Revert reference/ fixes

* Revert Changelog.txt typos

* Disable the SkyLakeX DGEMMITCOPY kernel as well

as a stopgap measure for numpy/numpy#13401 as mentioned in OpenMathLib#1955

* Disable DGEMMINCOPY as well for now

OpenMathLib#1955

* init

* Fix errors in cpu enumeration with glibc 2.6

for OpenMathLib#2114

* Change two http links to https

Closes OpenMathLib#2109

* remove redundant code OpenMathLib#2113

* Set up CI with Azure Pipelines

[skip ci]

* TST: add native POWER8 to CI

* add native POWER8 testing to
Travis CI matrix with ppc64le
os entry

* Update link to IBM MASS library, update cpu support status

* first try migrating one of the arm builds from travis

* fix tabbing in azure commands

* Update azure-pipelines.yml

take out offending lines (although stolen from https://github.com/conda-forge/opencv-feedstock azure-pipelines fiie)

* Update azure-pipelines.yml

* Update azure-pipelines.yml

* Update azure-pipelines.yml

* Update azure-pipelines.yml

* DOC: Add Azure CI status badge

* Add ARMV6 build to azure CI setup (OpenMathLib#2122)

using aytekinar's Alpine image and docker script from the Travis setup

[skip ci]

* TST: Azure manylinux1 & clean-up

* remove some of the steps & comments
from the original Azure yml template

* modify the trigger section to use
develop since OpenBLAS primarily uses
this branch; use the same batching
behavior as downstream projects NumPy/
SciPy

* remove Travis emulated ARMv6 gcc build
because this now happens in Azure

* use documented Ubuntu vmImage name for Azure
and add in a manylinux1 test run to the matrix

[skip appveyor]

* Add NO_AFFINITY to available options on Linux, and set it to ON

to match the gmake default. Fixes second part of OpenMathLib#2114

* Replace ISMIN and ISAMIN kernels on all x86_64 platforms (OpenMathLib#2125)

* Mark iamax_sse.S as unsuitable for MIN due to issue OpenMathLib#2116
* Use iamax.S rather than iamax_sse.S for ISMIN/ISAMIN on all x86_64 as workaround for OpenMathLib#2116

* Move ARMv8 gcc build from Travis to Azure

* Move ARMv8 gcc build from Travis to Azure

* Update .travis.yml

* Test drone CI

* install make

* remove sudo

* Install gcc

* Install perl

* Install gfortran and add a clang job

* gfortran->gcc-gfortran

* Switch to ubuntu and parallel jobs

* apt update

* Fix typo

* update yes

* no need of gcc in clang build

* Add a cmake build as well

* Add cmake builds and print options

* build without lapack on cmake

* parallel build

* See if ubuntu 19.04 fixes the ICE

* Remove qemu armv8 builds

* arm32 build

* Fix typo

* TST: add SkylakeX AVX512 CI test

* adapt the C-level reproducer code for some
recent SkylakeX AVX512 kernel issues, provided
by Isuru Fernando and modified by Martin Kroeker,
for usage in the utest suite

* add an Intel SDE SkylakeX emulation utest run to
the Azure CI matrix; a custom Docker build was required
because Ubuntu image provided by Azure does not support
AVX512VL instructions

* Add option USE_LOCKING for single-threaded build with locking support

for calling from concurrent threads

* Add option USE_LOCKING for single-threaded build with locking support

* Add option USE_LOCKING for SMP-like locking in USE_THREAD=0 builds

* Add option USE_LOCKING but keep default settings intact

* Remove unrelated change

* Do not try ancient PGI hacks with recent versions of that compiler

should fix OpenMathLib#2139

* Build and run utests in any case, they do their own checks for fortran availability

* Add softfp support in min/max kernels

fix for OpenMathLib#1912

* Revert "Add softfp support in min/max kernels"

* Separate implementations of AMAX and IAMAX on arm

As noted in OpenMathLib#1912 and comment on OpenMathLib#1942, the combined implementation happens to "do the right thing" on hardfp, but cannot return both value and index on softfp where they would have to share the return register

* Ensure correct output for DAMAX with softfp

* Use generic kernels for complex (I)AMAX to support softfp

* improved zgemm power9 based on power8

* upload thread safety test folder

* hook up c++ thread safety test (main Makefile)

*  add c++ thread test option to Makefile.rule

* Document NO_AVX512 

for OpenMathLib#2151

* sgemm pipeline improved, zgemm rewritten without inner packs, ABI lxvx v20 fixed with vs52

* Fix detection of AVX512 capable compilers in getarch

21eda8b introduced a check in getarch.c to test if the compiler is capable of
AVX512. This check currently fails, since the used __AVX2__ macro is only
defined if getarch itself was compiled with AVX2/AVX512 support. Make sure this
is the case by building getarch with -march=native on x86_64. It is only
supposed to run on the build host anyway.

* c_check: Unlink correct file

* power9 zgemm ztrmm optimized

* conflict resolve

* Add gfortran workaround for ABI violations in LAPACKE

for OpenMathLib#2154 (see gcc bug 90329)

* Add gfortran workaround for ABI violations

for OpenMathLib#2154 (see gcc bug 90329)

* Add gfortran workaround for potential ABI violation 

for OpenMathLib#2154

* Update fc.cmake

* Remove any inadvertent use of -march=native from DYNAMIC_ARCH builds

from OpenMathLib#2143, -march=native precludes use of more specific options like -march=skylake-avx512 in individual kernels, and defeats the purpose of dynamic arch anyway.

* Avoid unintentional activation of TLS code via USE_TLS=0

fixes OpenMathLib#2149

* Do not force gcc options on non-gcc compilers

fixes compile failure with pgi 18.10 as reported on OpenBLAS-users

* Update Makefile.x86_64

* Zero ecx with a mov instruction

PGI assembler does not like the initialization in the constraints.

* Fix mov syntax

* new sgemm 8x16

* Update dtrmm_kernel_16x4_power8.S

* PGI compiler does not like -march=native

* Fix build on FreeBSD/powerpc64.

Signed-off-by: Piotr Kubaj <pkubaj@anongoth.pl>

* Fix build for PPC970 on FreeBSD pt. 1

FreeBSD needs DCBT_ARG=0 as well.

* Fix build for PPC970 on FreeBSD pt.2

FreeBSD needs those macros too.

* cgemm/ctrmm power9

* Utest needs CBLAS but not necessarily FORTRAN

* Add mingw builds to Appveyor config

* Add getarch flags to disable AVX on x86

(and other small fixes to match Makefile behaviour)

* Make disabling DYNAMIC_ARCH on unsupported systems work

needs to be unset in the cache for the change to have any effect

* Mingw32 needs leading underscore on object names

(also copy BUNDERSCORE settings for FORTRAN from the corresponding Makefile)
@isuruf
Copy link
Contributor

isuruf commented Jul 19, 2019

Here's a C reproducer for the problem @andreasnoack mentioned. I used @fenrus75's code and modified it.

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <malloc.h>
#include <stdint.h>
#include <assert.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <f77blas.h>


double *m1, *m2;


#define PARAM 8


double *identity(int n)
{
	int i;
	double *m;
	m = calloc(n, n * sizeof(double));
	for (i = 0; i < n; i++)
		m[i + i * n] = 1.0;
	return m;
}

double * ones(int n)
{
	int i;
	double *m;
	m = calloc(n, n * sizeof(double));
	memset(m, 0, n*n*sizeof(double));
	for (i = 0; i < n*n; i++)
		m[i] = 1.0;
	return m;
}


void print_matrix(double *m, int n)
{
	int x,y;
	printf("----------------------------------------------\n");
	for (y = 0; y < n; y++) {
		for (x = 0; x < n; x++)
			printf("%5.2f\t", m[y * n + x]);
		printf("\n");
	}
	printf("==============================================\n");
}

void check_lower_left(double *m, int n)
{
	int x, y;
	for (y = 0; y < n ; y++) {
		for (x = 0; x < n; x++) {
			if (x < y && m[x + y * n] != 0.0)
				printf("Failure: value not 0.0 at %i, %i\n", y, x);
			if (x >= y && m[x + y * n] != 1.0)
				printf("Failure: value not 1.0 at %i, %i\n", y, x);
		}
	}
}



int main(int argc, char **argv)
{
	m1 = ones(PARAM);
	m2 = identity(PARAM);

	print_matrix(m1, PARAM);
	print_matrix(m2, PARAM);

    char SIDE='R';
    char UPLO='L';
    char TRANS='N';
    char DIAG='N';

    int p = PARAM;
    double alpha = 1.0;

	dtrmm_(&SIDE, &UPLO, &TRANS, &DIAG, &p, &p, &alpha, m1, &p, m2, &p);

	print_matrix(m2, PARAM);
	check_lower_left(m2, PARAM);
    free(m1);
    free(m2);
}

@isuruf
Copy link
Contributor

isuruf commented Jul 19, 2019

Looks like DTRMMKERNEL=dtrmm_kernel_4x8_haswell.c needs DGEMMI*COPY=gemm_*copy_4.c variants and DGEMMKERNEL=dgemm_kernel_4x8_skylakex.c needs DGEMMI*COPY=gemm_*copy_8.c variants.

@martin-frbg
Copy link
Collaborator

Ouch. That would certainly explain why each code in isolation "looks" correct. In retrospect, fenrus75's code started out as a 16x2 kernel with both DGEMM and DTRMM variants - all consistent, but unfortunately based on an unused and broken Haswell kernel from wernsaar's early work. Thereafter, fenrus75's work focused on the DGEMM kernel with no TRMM counterpart, and the current "4x8" must have already morphed into an 8x8 while keeping the original name at the time of the PR. I recall now that you already flagged the apparent discrepancy in a comment some weeks ago but I obviously did not understand the implications.

@isuruf
Copy link
Contributor

isuruf commented Jul 19, 2019

and the current "4x8" must have already morphed into an 8x8 while keeping the original name at the time of the PR

Thanks for confirming. I suspected this was the case, but the code here is way over my head to make such a claim. What's more confusing is that the value of DGEMM_DEFAULT_UNROLL_M is 4 in param.h for SkyLakeX x86_64.

@isuruf
Copy link
Contributor

isuruf commented Jul 19, 2019

I can try writing naive C code for a DTRMM 8x8 kernel and we can ask people to test. To do that I need to know what the difference between dtrmm_kernel and plain dtrmm_ is. They don't do the same thing right?

@martin-frbg
Copy link
Collaborator

martin-frbg commented Jul 20, 2019

Plain dtrmm_ is just the interface (in interface/trsm.c) that forwards to the appropriate function (either the gemm_thread_(m/n) subdivider for multithreaded, or directly to e.g. dtrmm_LNUU from driver/level3.
(BTW my comment should not be taken as confirmation of your theory, just that it looks plausible to me - I consider myself a sorcerer's apprentice at best. It is probably not a coincidence however that all other cpu-specific codes use identical MxN combinations for TRMM and GEMM.
Playing with DGEMM_DEFAULT_UNROLL_M or removing what "looked like" 8x8 code from the dgemm kernel has not led me to anything useful yet)

@martin-frbg
Copy link
Collaborator

Rechecking with BLAS-Tester - both DGEMM and DTRMM pass with fenrus75' code (as do TRSM and SYRK/SYR2K), only DSYMM fails - for any m greater than 7.

@isuruf
Copy link
Contributor

isuruf commented Jul 20, 2019

DTRMM fails with the code at #1955 (comment)

@martin-frbg
Copy link
Collaborator

martin-frbg commented Jul 20, 2019

Indeed it does. But I think I have hacked it into submission for now by "judicious" 😄 use of #if 0 whereever a loop looked related to the "greater than 7" observation. (Now to see - if I can - where the actual problem lies...)
dgemm_kernel_4x8_skylakex.txt

@isuruf
Copy link
Contributor

isuruf commented Jul 21, 2019

What did you use for DGEMMI*COPY, DTRMMKERNEL ?

@martin-frbg
Copy link
Collaborator

*copy_8_skylakex.c and the haswell dtrmm kernel - i.e. just the current KERNEL.SKYLAKEX with everything uncommented again as originally intended.

@martin-frbg
Copy link
Collaborator

Oh joy - I cannot reproduce any part of yesterday's "success" - in fact that mangled version fails just about any BLAS3 test. I can only assume I managed to slip myself a non-AVX512 build somehow.
Guess I should just release 0.3.7 as it is today and take a break from the project (apart from merging PRs obviously)

@martin-frbg
Copy link
Collaborator

This combo does appear to work -
dgemm_kernel_4x8_skylakex.txt
KERNEL.SKYLAKEX.txt
but I suspect the performance of the brutally stripped-down dgemm kernel may actually be below the Haswell one now. (Cannot check this as I am using SDE on older hardware)

@martin-frbg
Copy link
Collaborator

Not surprisingly, that hacked version performs markedly worse in the DGEMM benchmark than the regular Haswell microkernel, now that I have the hardware to test this.

@martin-frbg martin-frbg modified the milestones: 0.3.7, 0.3.8 Aug 11, 2019
@martin-frbg
Copy link
Collaborator

I have failed to come up with a proper solution to this so far, so will release 0.3.7 in its current state i.e. with the AVX512 DGEMM kernel parts disabled.

@martin-frbg
Copy link
Collaborator

Fixed now by wjc404's new AVX512 DGEMM kernel from #2286

@isuruf
Copy link
Contributor

isuruf commented Oct 24, 2019

@staticfloat, @matthew-brett, can you check that the develop branch works?

@matthew-brett
Copy link
Contributor

@tylerjereddy - sorry - I'm a bit out of the loop on Numpy CI - are you the right person to ask about this?

@tylerjereddy
Copy link
Contributor

I'll try to bump OpenBLAS in the MacPython ecosystem soon so we can see if we can leverage the new AVX512 kernel(s) that were previously disabled. Note that both NumPy and SciPy decided not to guard against this in CI, since the OpenBLAS CI currently includes an appropriate Intel Emulator regression test for the matter.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests