Openblas: DGEMM regression on SkylakeX
It looks like https://github.com/xianyi/OpenBLAS/commit/45fe8cb0c5d06f890913e86078cb48ac379c65dc 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 544b069e85254d41699afde16e2e81c123cb5f28) 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.
staticfloat
All 85 comments
It would be good to have a point release with a fix for this soon, if possible.
ViralBShah
on 9 Jan 2019
We're happy to test a patch on our entire testsuite and report before the release is tagged, in case other issues are uncovered.
ViralBShah
on 9 Jan 2019
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."
brianborchers
on 9 Jan 2019
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.
ViralBShah
on 9 Jan 2019
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.
brianborchers
on 9 Jan 2019
There are lots of other tests failing, it turns out. I am trying to find a simple testcase to reproduce.
ViralBShah
on 9 Jan 2019
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.
staticfloat
on 9 Jan 2019
Original issue was #1643 , too bad this did not get caught in 0.3.4 but I still do not have the hardware.
martin-frbg
on 9 Jan 2019
@staticfloat is it possible to drill problem down to copy vs kernel routines?
brada4
on 9 Jan 2019
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.
martin-frbg
on 9 Jan 2019
@keno can we give @martin-frbg access to our Skylake machine - Antarctic?
ViralBShah
on 9 Jan 2019
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.
staticfloat
on 9 Jan 2019
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.
martin-frbg
on 9 Jan 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
You could copy respective options from known to work KERNEL.HASWEL
brada4
on 9 Jan 2019
@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
on 9 Jan 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
So DGEMMKERNEL is definitely broken and needs to be commented out ? Will probably need help from @fenrus75 who wrote this kernel.
martin-frbg
on 10 Jan 2019
back to test case.. is there anything I can easily test this with?
fenrus75
on 10 Jan 2019
(also @martin-frbg send me an email and we'll talk about skx hardware)
fenrus75
on 10 Jan 2019
ping? is there any way to reproduce a failure ?
fenrus75
on 24 Jan 2019
@fenrus75 do you need access to a SkylakeX machine? I can get you SSH access on a Linux SkylakeX machine if that would help.
staticfloat
on 24 Jan 2019
I work for Intel... in my cube alone I have 3 skylakeX's ;-)
I just need code that shows it go wrong so that I can debug and fix
On Thu, Jan 24, 2019 at 1:43 PM Elliot Saba notifications@github.com
wrote:
@fenrus75 https://github.com/fenrus75 do you need access to a SkylakeX
machine? I can get you SSH access on a Linux SkylakeX machine if that would
help.—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/xianyi/OpenBLAS/issues/1955#issuecomment-457368306,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ABPeFf1n8wk3zqRJSIL5VMfdpnUzWEBVks5vGikagaJpZM4Z2an-
.
fenrus75
on 24 Jan 2019
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.
staticfloat
on 25 Jan 2019
It's actually relatively complicated. The Julia code is calling a series of LAPACK calls to find the pseudoinverse of a matrix and then testing that the result satisfies the properties of the pseudoinverse. That test of the properties of a pseudoinverse is failing. This could be because the pseudoinverse is computed incorrectly (due to incorrect results from a BLAS call deep inside LAPACK) or because of incorrect computations in checking the pseudoinverse properties.
Unfortunately, the particular matrix used in the Julia tests is very badly conditioned, so there isn't really a single "correct" answer- lots and lots of floating point matrices come sufficiently close to satisfying these pseudoinverse properties and could be considered "correct" answers. For example, the might (1,1) entry in the pseudoinverse matrix might very well be -1.0e6 or +1e6 in two different matrices that satisfy the pseudoinverse conditions. This can easily happen due to changing the order of FMA's and the non-associativity of floating point arithmetic.
I would start by making sure that the Skylake build of the current working version correctly passes its own build-time tests. If any of those fails, I'd start to debug from that rather than the Julia test. If not, I'd then look at the BLAS/LAPACK routines used in the Julia computations one by one, checking that they produce correct answers or similarly sized but well-conditioned test problems with well-known answers.
brianborchers
on 25 Jan 2019
It's likely the Julia is using the LAPACK routine DGESDD to compute the SVD. However, because this matrix is so badly conditioned, most of the singular values computed by the routine will be essentially random numbers of about 1.0e-16 times the size of the largest singular value. You could use instead a well-conditioned matrix for which the singular values should come out consistently.
brianborchers
on 25 Jan 2019
if I know the exact arguments to the blas function I might be able to make a reproducer that's simpler
my own tests run many different sizes, but with pretty conventional arguments for all the non-size arguments to the gemm functions
fenrus75
on 25 Jan 2019
Does DGEMM work OK for a variety of matrix sizes with the Skylake build? If so, I'd tend to believe that DGEMM was fine but that some other BLAS routine was broken. I think it's unlikely that any of the LAPACK routines are broken. DGESDD calls a bunch of LAPACK auxilliary routines, and they in turn can call a bunch of BLAS routines.
Have you tried testing DGESDD on some random but reasonably well-conditioned test cases? For such matrices, you could compare the DGESDD computed singular values with singular values computed using a known good BLAS and compare them to test whether DGESDD was working correctly.
Again, the reason that you can't do this with the very badly conditioned matrix in the Julia test is that the problem is so badly conditioned that absolutely any difference in the order of arithmetic operations can result in vastly different answers (and you probably won't get reproducible answers from run to run either.)
The way that Julia's test handles this ill-conditioned matrix in the test is to confirm that the returned pseudoinverse satisfies the Moore-Penrose properties, not that it matches a particular known answer. Testing those Moore-Penrose properties requires more calls to BLAS, and it could be one of these that is where the bug appears.
brianborchers
on 25 Jan 2019
I'll run even more tests... but dgemm takes many different options in terms of ordering of data, transposing, strides, beta, etc etc and there can be subtle bugs in that.
fenrus75
on 25 Jan 2019
Looking at the code for DGESDD, in the calls to DGEMM, none of the matrices are transposed, all of the alpha values are 1, and all of the beta values are 0. The M, N, K size parameters do vary. There aren't any strides in DGEMM.
brianborchers
on 25 Jan 2019
It might be easier to find a simple reproducible test case in one of the many other Julia linear algebra tests that is failing.
brianborchers
on 25 Jan 2019
https://github.com/Reference-LAPACK/lapack/issues/316 is probably unrelated (except for DGESDD with an ill-conditioned matrix) but apparently netlib switches to a less stable divide&conquer algorithm at a certain matrix size. @fenrus75 I assume brianborchers' test case from #1643 still passes with your latest changes ?
martin-frbg
on 25 Jan 2019
I opened issue #316 for LAPACK the other day and I gave further details on about the how the errors in the singular values from GESDD can be arbitrarily bad in my corresponding Octave bug report. This is in line with brianborchers's comment above. That said, GESVD seems to be much better at computing the singular values for such difficult matrices matrices. An easy way to see this is to shift the problematic matrix to one of its own eigenvalues, so that the result is singular. For the example I tested, the maximum singular value was roughly 6e+25 (so huge). However, GESVD still properly returned that the smallest singular value was roughly 1e-23 whereas GESDD returned roughly 1e+9 as the smallest singular value!
At the moment, we've yet to receive confirmation from the LAPACK people whether the specific behavior I noticed is due to a bug in the code or a simply a limitation of the divide and conquer algorithm implemented in GESDD. In any case, it seems avoiding GESDD in favor of GESVD would be prudent, at least for the time being.
t-mitchell
on 25 Jan 2019
This morning, I've spent some time on the original pinv issue mentioned in the top post and my partial conclusion so far is that it's different from https://github.com/Reference-LAPACK/lapack/issues/316. In our example here, the singular values are roughly the same. The differences among the very small values don't matter since they are ignored in the pinv computation. Instead, the problem here seems to be the singular vectors. When computed with dgesdd, the U matrix is not even close to being orthogonal. Interestingly, it doesn't seem to be an issue with the bidiagonal DnC solver dbdsdc. It produces fine results based on the bidiagonal matrix from the bidiagonal reduction of the Hilbert matrix.
There are a lot of different branches in LAPACKs SVD solvers depending on the aspect ratio of the matrix, the available workspace, and how many of the singular vectors are requested so it's a bit complicated to figure out exactly where the DnC and QR based solvers differ but I think the issue here is related to how the elementary (Householder) reflectors are applied. Try computing just the QR of the Hilbert matrix and verify that Q is not even close to orthogonal.
julia> Matrix(Q)'*Matrix(Q)
100×100 Array{Float64,2}:
2.72358 0.836204 0.541382 1.56653 0.723944 -0.631045 0.393747 0.415072 … -0.379568 -0.389309 -0.171272 0.00764817 -0.00274127 -0.00142251 -0.00329317
0.836204 5.22922 -0.199689 1.49522 -0.84534 0.986203 0.12605 0.806151 -0.0323286 0.254786 0.123991 0.0630916 0.00275348 -0.0165985 0.0317004
0.541382 -0.199689 4.48407 0.626691 -0.34602 -1.01427 1.26365 1.36967 0.115609 -1.55789 -0.329611 -0.0178857 0.00210697 0.0052537 -0.0105639
1.56653 1.49522 0.626691 7.09471 -0.427999 -1.15392 -0.33465 1.88252 -0.620717 -0.499573 -0.281308 0.0625795 0.000947451 -0.0187098 0.0225606
0.723944 -0.84534 -0.34602 -0.427999 2.06863 -0.262304 -0.646683 -0.647517 -0.412086 0.0792543 0.170981 -0.000300248 0.000309771 3.87078e-5 5.30124e-5
-0.631045 0.986203 -1.01427 -1.15392 -0.262304 2.62505 0.530753 -0.292458 … -0.473683 0.820261 0.0993953 0.00900056 -0.00447621 -0.00160283 0.00109047
0.393747 0.12605 1.26365 -0.33465 -0.646683 0.530753 3.53951 0.709963 -0.122997 -0.385725 -0.475491 -0.00719675 -0.00788346 0.00499399 -0.0110819
0.415072 0.806151 1.36967 1.88252 -0.647517 -0.292458 0.709963 3.01312 0.379387 0.314183 -0.238069 0.0652777 -0.0218099 -0.0138893 0.0058541
-0.307765 0.0203183 -0.211825 -1.00913 0.0379945 0.690239 -0.0336199 -0.210315 -0.132379 0.291022 0.0796898 -2.73219e-17 7.04731e-19 -6.07153e-18 -3.44234e-18
4.47614e-5 0.196967 -0.42677 -0.326762 0.274565 0.46265 -0.323948 -0.17351 -0.359841 0.317091 0.114012 -9.1073e-18 -4.87891e-18 -3.46945e-18 -1.33357e-17
0.0258762 0.441137 -0.530672 0.339885 -0.399229 0.210467 0.341043 0.0941537 … -0.275301 0.850859 0.0709662 5.63785e-18 7.11508e-18 3.68629e-18 -1.33357e-17
⋮ ⋮ ⋱ ⋮
-0.258421 0.71181 -0.861926 -0.847345 -0.164887 0.175783 -0.523864 -0.665899 0.395669 0.962552 0.25602 -0.151686 0.0395261 0.0460646 -0.0873332
-0.0823063 1.42406 -1.63907 -0.25803 -0.154554 0.396869 -0.318054 -0.401972 … -0.0260761 -0.0558141 0.993167 0.160552 -0.261107 0.138922 -0.0218541
-0.191812 0.594102 -0.334123 -0.327639 -0.0986163 0.0699678 0.162216 -0.150285 -0.24146 1.11374 -0.927441 -0.043752 0.154542 -0.0486787 0.070782
-0.239292 -0.283033 0.450502 0.207118 -0.21345 -0.42476 0.153402 0.0297775 0.691516 -0.884044 0.281347 0.0609737 -0.00230397 -0.061953 0.0510808
-0.379568 -0.0323286 0.115609 -0.620717 -0.412086 -0.473683 -0.122997 0.379387 6.04665 -0.280263 0.691108 0.203946 -0.0852124 -0.0542395 -0.0966923
-0.389309 0.254786 -1.55789 -0.499573 0.0792543 0.820261 -0.385725 0.314183 -0.280263 9.26107 0.503453 -0.122143 -0.15845 0.2411 0.0188728
-0.171272 0.123991 -0.329611 -0.281308 0.170981 0.0993953 -0.475491 -0.238069 … 0.691108 0.503453 5.03552 0.120413 -0.034236 0.0439063 -0.198219
0.00764817 0.0630916 -0.0178857 0.0625795 -0.000300248 0.00900056 -0.00719675 0.0652777 0.203946 -0.122143 0.120413 1.0 5.63785e-17 6.59195e-17 -4.16334e-17
-0.00274127 0.00275348 0.00210697 0.000947451 0.000309771 -0.00447621 -0.00788346 -0.0218099 -0.0852124 -0.15845 -0.034236 5.63785e-17 1.0 5.20417e-18 1.04083e-16
-0.00142251 -0.0165985 0.0052537 -0.0187098 3.87078e-5 -0.00160283 0.00499399 -0.0138893 -0.0542395 0.2411 0.0439063 6.59195e-17 5.20417e-18 1.0 -6.50521e-17
-0.00329317 0.0317004 -0.0105639 0.0225606 5.30124e-5 0.00109047 -0.0110819 0.0058541 -0.0966923 0.0188728 -0.198219 -4.16334e-17 1.04083e-16 -6.50521e-17 1.0
I'll continue the investigation.
andreasnoack
on 25 Jan 2019
Just a short update here in case other people are also trying to figure out what is going on. I'm pretty sure the issue is in dormqr, i.e. the blocked routine to multiply with a matrix of elementary (Householder) reflectors. Given the same reflectors, I get an orthogonal matrix if I use the unblocked version dorm2r instead. This also fits well with the difference between the QR and DnC based solvers since the QR based solver dgesvd only calls dormqr if the work space isn't sufficiently large. This branch isn't taken since we query for the optimal work space size. In contrast, the DnC based solver dgesdd calls dormqr even when using the optimal work space size (as I read the source) so this could explain the difference.
I'm still not sure exactly where things go wrong inside dormqr. It again calls dlarft and dlarfb and the latter is the only one that calls dgemm so that is currently my candidate.
andreasnoack
on 25 Jan 2019
BINGO. The issue is dtrmm when multiplying from the right
julia> n = 7;
julia> BLAS.trmm('R', 'L', 'N', 'N', 1.0, ones(n, n), Matrix{Float64}(I, n, n))
7×7 Array{Float64,2}:
1.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 1.0 0.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0
julia> n = 8;
julia> BLAS.trmm('R', 'L', 'N', 'N', 1.0, ones(n, n), Matrix{Float64}(I, n, n))
8×8 Array{Float64,2}:
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0
1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
Isn't this case supposed to be covered by the OpenBLAS tests?
andreasnoack
on 25 Jan 2019
Well apparently not ... could be due to an unlucky choice of N values for the lapack dblat3 test (... 7, 31, 63 in dblat3.dat) unless it got overlooked by those with the hardware (the CI tests use Nehalem kernels)
martin-frbg
on 25 Jan 2019
I'll build a C based testcase and figure out from there.... going to be an
interesting weekend project
On Fri, Jan 25, 2019 at 7:18 AM Martin Kroeker notifications@github.com
wrote:
Well apparently not ... could be due to an unlucky choice of N values for
the lapack dblat3 test (... 7, 31, 63 in dblat3.dat) unless it got
overlooked by those with the hardware (the CI tests use Nehalem kernels)—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/xianyi/OpenBLAS/issues/1955#issuecomment-457606054,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ABPeFQvfDuy4suiegA-LfDTt-VChKMOFks5vGyA5gaJpZM4Z2an-
.
fenrus75
on 25 Jan 2019
#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 <cblas.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] != 1.0)
printf("Failure: value not 1.0 at %i, %i\n", x,y);
if (x > y && m[x + y * n] != 0.0)
printf("Failure: value not 0.0 at %i, %i\n", x,y);
}
}
}
int main(int argc, char **argv)
{
m1 = ones(PARAM);
m2 = identity(PARAM);
print_matrix(m1, PARAM);
print_matrix(m2, PARAM);
cblas_dtrmm(CblasRowMajor, CblasRight, CblasLower, CblasNoTrans, CblasNonUnit, PARAM, PARAM, 1.0, m1, PARAM, m2, PARAM);
print_matrix(m2, PARAM);
check_lower_left(m2, PARAM);
}
basic start of a test case in previous comment; have not succeeded in making it fail yet
fenrus75
on 26 Jan 2019
Hm. Something more subtle is going on here. I'm now unable to reproduce the minimal example but I still get that Matrix(qr(randn(n,n)).Q) |> t -> t't is the identity if n<=56 and completely off for n>56. I'll take another look.
andreasnoack
on 26 Jan 2019
it is 57 or 64 that breaks?
(does matter since there's codepaths that are different between those)
fenrus75
on 27 Jan 2019
Right now it fails consistently at 57 in Julia but not in Fortran.
andreasnoack
on 27 Jan 2019
I assume Julia uses C matrix ordering not Fortran.. in which case the issue
is most likely on the incopy side somehow
On Sun, Jan 27, 2019 at 12:16 AM Andreas Noack notifications@github.com
wrote:
Right now it fails consistently at 57 in Julia but not in Fortran.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/xianyi/OpenBLAS/issues/1955#issuecomment-457898484,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ABPeFWMTzrvPY8HfMy1LrgQ_8rNvUL4nks5vHWBogaJpZM4Z2an-
.
fenrus75
on 27 Jan 2019
Glancing through the ncopy file while already half asleep, is this correct:
https://github.com/xianyi/OpenBLAS/blob/89b60dab8ad21a0cc6320cbd9fcd603c4c4bfc81/kernel/x86_64/dgemm_ncopy_8_skylakex.c#L377-L380
(the assigning from ctemp01..._3_..._2_...4) ?
martin-frbg
on 27 Jan 2019
ohhhh no it does not appear to be!
note that this is the same as in the generic/ code that this function
originally came from
(
https://github.com/xianyi/OpenBLAS/blob/89b60dab8ad21a0cc6320cbd9fcd603c4c4bfc81/kernel/generic/gemm_ncopy_8.c#L379
)
On Sun, Jan 27, 2019 at 2:40 PM Martin Kroeker notifications@github.com
wrote:
Glancing through the ncopy file while already half asleep, is this correct:
https://github.com/xianyi/OpenBLAS/blob/89b60dab8ad21a0cc6320cbd9fcd603c4c4bfc81/kernel/x86_64/dgemm_ncopy_8_skylakex.c#L377-L380
(the assigning from ctemp01...3...2...4) ?—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/xianyi/OpenBLAS/issues/1955#issuecomment-457961190,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ABPeFWiFsnPKW59A9SXOla7q2M3_wGvSks5vHirQgaJpZM4Z2an-
.
fenrus75
on 28 Jan 2019
hmm it's not so obvious looking at it more; the code higher up in the 8x unroll case also has equivalent interleaving
fenrus75
on 28 Jan 2019
I assume Julia uses C matrix ordering not Fortran
Julia has column major arrays like Fortran and we also use the Fortran API when calling BLAS.
andreasnoack
on 28 Jan 2019
Any progress on this ? I assume the DGEMM problem is not reproducible with something comparatively simple like the xdl3blastst from BLAS-Tester (or ATLAS originally) ?
martin-frbg
on 20 Feb 2019
bin/xdl3blastst
--------------------------------- GEMM ----------------------------------
TST# A B M N K ALPHA LDA LDB BETA LDC TIME MFLOP SpUp TEST
==== = = ==== ==== ==== ===== ==== ==== ===== ==== ===== ===== ==== =====
0 N N 100 100 100 1.0 1000 1000 1.0 1000 0.17 11.5 1.00 -----
0 N N 100 100 100 1.0 1000 1000 1.0 1000 0.01 196.0 17.10 PASS
1 N N 200 200 200 1.0 1000 1000 1.0 1000 0.00 5753.0 1.00 -----
1 N N 200 200 200 1.0 1000 1000 1.0 1000 0.00 38836.1 6.75 PASS
2 N N 300 300 300 1.0 1000 1000 1.0 1000 0.01 5031.7 1.00 -----
2 N N 300 300 300 1.0 1000 1000 1.0 1000 0.00 143622.3 28.54 PASS
3 N N 400 400 400 1.0 1000 1000 1.0 1000 0.04 3657.1 1.00 -----
3 N N 400 400 400 1.0 1000 1000 1.0 1000 0.00 171853.7 46.99 PASS
4 N N 500 500 500 1.0 1000 1000 1.0 1000 0.06 4431.6 1.00 -----
4 N N 500 500 500 1.0 1000 1000 1.0 1000 0.00 326862.8 73.76 PASS
5 N N 600 600 600 1.0 1000 1000 1.0 1000 0.09 4957.1 1.00 -----
5 N N 600 600 600 1.0 1000 1000 1.0 1000 0.00 415868.6 83.89 PASS
6 N N 700 700 700 1.0 1000 1000 1.0 1000 0.14 5070.3 1.00 -----
6 N N 700 700 700 1.0 1000 1000 1.0 1000 0.00 379189.8 74.79 PASS
7 N N 800 800 800 1.0 1000 1000 1.0 1000 0.20 5162.6 1.00 -----
7 N N 800 800 800 1.0 1000 1000 1.0 1000 0.00 412422.4 79.89 PASS
8 N N 900 900 900 1.0 1000 1000 1.0 1000 0.26 5556.7 1.00 -----
8 N N 900 900 900 1.0 1000 1000 1.0 1000 0.00 528091.1 95.04 PASS
9 N N 1000 1000 1000 1.0 1000 1000 1.0 1000 0.35 5749.0 1.00 -----
9 N N 1000 1000 1000 1.0 1000 1000 1.0 1000 0.00 463946.0 80.70 PASS
On Wed, Feb 20, 2019 at 5:35 AM Martin Kroeker notifications@github.com
wrote:
Any progress on this ? I assume the DGEMM problem is not reproducible with
something comparatively simple like the xdl3blastst from BLAS-Tester (or
ATLAS originally) ?—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/xianyi/OpenBLAS/issues/1955#issuecomment-465575709,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ABPeFcGPeXlCLjpXI4pUnlIjC92c5ypiks5vPU8QgaJpZM4Z2an-
.
fenrus75
on 20 Feb 2019
Thanks. Same for TRMM (xdl3blastst -R trmm) and varied values of M,N,K I guess... 90x speedup compared to netlib looks nice though
martin-frbg
on 20 Feb 2019
I haven't been able to reproduce this in Fortran. Any ideas to why this could be? Memory alignment differences?
andreasnoack
on 20 Feb 2019
Perhaps running some of the recent test cases under valgrind would tell if there are any unintended accesses ? (Would have thought something like this much more likely with the old hand-crafted assembly rather than those newfangled intrinsics though...)
martin-frbg
on 20 Feb 2019
I don't know whether this would help to reproduce the problem, but I think we just hit the same problem in Numpy: https://github.com/numpy/numpy/issues/13401
The result of the SVD is way off, and reliably so, it seems, at the Numpy level, when using the SkylakeX kernel.
matthew-brett
on 29 Apr 2019
Here's a C reproducer for the issue @matthew-brett mentioned. It fails with openblas 0.3.6 and SkyLakeX. (Doesn't print 0.0)
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <lapacke.h>
#define SIZE 100
double s[SIZE];
double u[SIZE*SIZE];
double vt[SIZE*SIZE];
double X[SIZE*SIZE];
double superb[SIZE];
double tmp[SIZE*SIZE];
double m[SIZE*SIZE];
int main(int argc, char **argv)
{
srand(0);
for (int i = 0; i < SIZE*SIZE; i++) {
m[i] = (rand()+0.0)/RAND_MAX * 10;
tmp[i] = m[i];
}
int info = LAPACKE_dgesvd( LAPACK_ROW_MAJOR, 'A', 'A', SIZE, SIZE, m, SIZE,
s, u, SIZE, vt, SIZE, superb);
for (int i = 0; i < SIZE; i++) {
for (int j = 0; j < SIZE; j++) {
u[i*SIZE+j] = u[i*SIZE+j]*s[j];
}
}
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
SIZE, SIZE, SIZE, 1, u, SIZE, vt, SIZE, 0, X, SIZE);
for (int i = 0; i < SIZE*SIZE; i++) {
X[i] = X[i] - tmp[i];
}
double norm = cblas_dnrm2(SIZE*SIZE, X, 1);
printf("%lf\n", norm);
}
cc @tylerjereddy
isuruf
on 4 May 2019
@isuruf thanks for the reproducer. If SDE is to be trusted, the culprit in that particular case is DGEMMITCOPY a.k.a. dgemm_tcopy_8_skylakex.c
martin-frbg
on 4 May 2019
Hmmm. The reproducer for numpy13401 actually passes with the AVX512 DGEMM kernel and ITCOPY both enabled (as does the xdl3blastst GEMM test from BLAS-Tester, but not the corresponding SYMM test as seen in #2029).
Reverting only the DGEMM kernel "fixed" SYMM but broke GEMM (which I should have noticed back then) and numpy13401. Disabling ITCOPY in addition makes 13401 and SYMM work but leaves the GEMM test still unhappy. Removing INCOPY as well, as suggested by staticfloat above makes all known failure cases work, so my fault for trying to go for a minimal change. :frowning_face:
Unfortunately I believe we are still no closer to finding the actual problem with the new AVX512 code.
@fenrus75 any suggestions ?
martin-frbg
on 5 May 2019
@isuruf I've been working a little on writing in a regression test based on your C reproducer (for i.e., emulation in CI here)--I have trouble with the #include <lapacke.h> line when I put code similar to yours in a utest/ C file for example.
I don't know much about the OpenBLAS test machinery though--that header should be available from building OpenBLAS with default make in time to run the utests, right? Or maybe not? Maybe utests is not the right spot for it?
tylerjereddy
on 11 May 2019
Lapacke and its header is built with C compiler on top of Fortran-compiled LAPACK libraries
brada4
on 11 May 2019
At the time the utests are run, lapacke.h is still buried in lapack-netlib/LAPACKE/include
martin-frbg
on 11 May 2019
Hmm, if I use that include path it helps, but I still have some issues. The new test code in my fork is visible here with includes.
I get a bunch of output like this, with the new test temporarily living in test_axpy.c:
In file included from ../lapack-netlib/LAPACKE/include/lapacke.h:151:0,
from test_axpy.c:38:
../lapack-netlib/LAPACKE/include/lapacke.h:12613:37: error: conflicting types for ‘sgetrf_’
#define LAPACK_sgetrf LAPACK_GLOBAL(sgetrf,SGETRF)
The definition of LAPACK_GLOBAL comes from lapacke_mangling.h, which lapacke.h tries to include - however, before installation this file is named lapacke_mangling_with_flags.h (not our idea, as this part of the code is imported from "netlib" Reference-LAPACK). Probably easiest to just copy that declaration into your test code.
martin-frbg
on 11 May 2019
This one seems to work for me (have not investigated why it does not hit 0.0 within DOUBLE_EPS on Haswell, but could be due to unrelated changes in my local tree)
test_axpy.txt
martin-frbg
on 11 May 2019
Thanks! I've been playing around with including the mangling header directly but didn't make much ground.
tylerjereddy
on 11 May 2019
So, the idea is to provide the function prototype in the absence of the header and then when it comes to test "runtime" the LAPACK(E) stuff will be built / all ready?
tylerjereddy
on 11 May 2019
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).
martin-frbg
on 12 May 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);
}
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.
isuruf
on 19 Jul 2019
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.
martin-frbg
on 19 Jul 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
on 19 Jul 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?
isuruf
on 19 Jul 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
on 20 Jul 2019
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.
martin-frbg
on 20 Jul 2019
DTRMM fails with the code at https://github.com/xianyi/OpenBLAS/issues/1955#issuecomment-513081029
isuruf
on 20 Jul 2019
Indeed it does. But I think I have hacked it into submission for now by "judicious" :smile: 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
martin-frbg
on 20 Jul 2019
What did you use for DGEMMI*COPY, DTRMMKERNEL ?
isuruf
on 21 Jul 2019
*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
on 21 Jul 2019
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
on 21 Jul 2019
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
on 21 Jul 2019
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
on 7 Aug 2019
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
on 11 Aug 2019
Fixed now by wjc404's new AVX512 DGEMM kernel from #2286
martin-frbg
on 20 Oct 2019
@staticfloat, @matthew-brett, can you check that the develop branch works?
isuruf
on 24 Oct 2019
@tylerjereddy - sorry - I'm a bit out of the loop on Numpy CI - are you the right person to ask about this?
matthew-brett
on 25 Oct 2019
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.
tylerjereddy
on 29 Nov 2019
Related issues
marker68
·
8Comments
nh2
·
9Comments
cswluo
·
5Comments
hetong007
·
5Comments
jakirkham
·
16Comments