Wrong results for small DTRMV on x86-64

[grisuthedragon]

Usually I use the POWER8 platform but today I found a bad error in the x86-64 code for the DTRMV routine.

I want to compute x <- A*x, where

• A is nonunit upper triangular,
• N = 1...64,
• LDA = 64,
• and INCX=1

and get wrong results on a Haswell (16 cores, Xeon CPU E5-2640 v3) based system using gcc (Ubuntu 5.4.0-6ubuntu1~16.04.5) if OpenBLAS is compiled via

make USE_OPENMP=1 

I compared the results with the Netlib implementation and obtained:

Correct RESULT: DTRMV(U,N,N,    1, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    2, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    3, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    4, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    5, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    6, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    7, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    8, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,    9, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   10, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   11, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   12, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   13, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   14, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   15, A,   64, X,     1)  MAXERR = .000000000000000D+00
Correct RESULT: DTRMV(U,N,N,   16, A,   64, X,     1)  MAXERR = .000000000000000D+00
Wrong RESULT  : DTRMV(U,N,N,   17, A,   64, X,     1)  MAXERR = .568434188608080D-13
Wrong RESULT  : DTRMV(U,N,N,   18, A,   64, X,     1)  MAXERR = .113686837721616D-12
Wrong RESULT  : DTRMV(U,N,N,   19, A,   64, X,     1)  MAXERR = .113686837721616D-12
Wrong RESULT  : DTRMV(U,N,N,   20, A,   64, X,     1)  MAXERR = .227373675443232D-12
Wrong RESULT  : DTRMV(U,N,N,   21, A,   64, X,     1)  MAXERR = .227373675443232D-12
Wrong RESULT  : DTRMV(U,N,N,   22, A,   64, X,     1)  MAXERR = .454747350886464D-12
Wrong RESULT  : DTRMV(U,N,N,   23, A,   64, X,     1)  MAXERR = .909494701772928D-12
Wrong RESULT  : DTRMV(U,N,N,   24, A,   64, X,     1)  MAXERR = .909494701772928D-12
Wrong RESULT  : DTRMV(U,N,N,   25, A,   64, X,     1)  MAXERR = .181898940354586D-11
Wrong RESULT  : DTRMV(U,N,N,   26, A,   64, X,     1)  MAXERR = .363797880709171D-11
Wrong RESULT  : DTRMV(U,N,N,   27, A,   64, X,     1)  MAXERR = .727595761418343D-11
Wrong RESULT  : DTRMV(U,N,N,   28, A,   64, X,     1)  MAXERR = .363797880709171D-11
Wrong RESULT  : DTRMV(U,N,N,   29, A,   64, X,     1)  MAXERR = .291038304567337D-10
Wrong RESULT  : DTRMV(U,N,N,   30, A,   64, X,     1)  MAXERR = .582076609134674D-10
Wrong RESULT  : DTRMV(U,N,N,   31, A,   64, X,     1)  MAXERR = .145519152283669D-10
Wrong RESULT  : DTRMV(U,N,N,   32, A,   64, X,     1)  MAXERR = .582076609134674D-10
Wrong RESULT  : DTRMV(U,N,N,   33, A,   64, X,     1)  MAXERR = .596598620663450D+06
Wrong RESULT  : DTRMV(U,N,N,   34, A,   64, X,     1)  MAXERR = .605368209769108D+06
Wrong RESULT  : DTRMV(U,N,N,   35, A,   64, X,     1)  MAXERR = .114251427897636D+07
Wrong RESULT  : DTRMV(U,N,N,   36, A,   64, X,     1)  MAXERR = .176011411015678D+07
Wrong RESULT  : DTRMV(U,N,N,   37, A,   64, X,     1)  MAXERR = .273520848158376D+07
Wrong RESULT  : DTRMV(U,N,N,   38, A,   64, X,     1)  MAXERR = .400535392194670D+07
Wrong RESULT  : DTRMV(U,N,N,   39, A,   64, X,     1)  MAXERR = .583437738276019D+07
Wrong RESULT  : DTRMV(U,N,N,   40, A,   64, X,     1)  MAXERR = .850833665131698D+07
Wrong RESULT  : DTRMV(U,N,N,   41, A,   64, X,     1)  MAXERR = .123493561359394D+08
Wrong RESULT  : DTRMV(U,N,N,   42, A,   64, X,     1)  MAXERR = .173391397036582D+08
Wrong RESULT  : DTRMV(U,N,N,   43, A,   64, X,     1)  MAXERR = .260484477387416D+08
Wrong RESULT  : DTRMV(U,N,N,   44, A,   64, X,     1)  MAXERR = .376086864422978D+08
Wrong RESULT  : DTRMV(U,N,N,   45, A,   64, X,     1)  MAXERR = .562860046618746D+08
Wrong RESULT  : DTRMV(U,N,N,   46, A,   64, X,     1)  MAXERR = .722429815166366D+08
Wrong RESULT  : DTRMV(U,N,N,   47, A,   64, X,     1)  MAXERR = .129301313570969D+09
Wrong RESULT  : DTRMV(U,N,N,   48, A,   64, X,     1)  MAXERR = .193137989170661D+09
Wrong RESULT  : DTRMV(U,N,N,   49, A,   64, X,     1)  MAXERR = .279222204411076D+09
Wrong RESULT  : DTRMV(U,N,N,   50, A,   64, X,     1)  MAXERR = .678662007793433D+09
Wrong RESULT  : DTRMV(U,N,N,   51, A,   64, X,     1)  MAXERR = .121269930403095D+10
Wrong RESULT  : DTRMV(U,N,N,   52, A,   64, X,     1)  MAXERR = .106155723183043D+10
Wrong RESULT  : DTRMV(U,N,N,   53, A,   64, X,     1)  MAXERR = .306408230139373D+10
Wrong RESULT  : DTRMV(U,N,N,   54, A,   64, X,     1)  MAXERR = .470622143023516D+10
Wrong RESULT  : DTRMV(U,N,N,   55, A,   64, X,     1)  MAXERR = .716298126556327D+10
Wrong RESULT  : DTRMV(U,N,N,   56, A,   64, X,     1)  MAXERR = .103315607770989D+11
Wrong RESULT  : DTRMV(U,N,N,   57, A,   64, X,     1)  MAXERR = .837779499609675D+10
Wrong RESULT  : DTRMV(U,N,N,   58, A,   64, X,     1)  MAXERR = .255092826195343D+11
Wrong RESULT  : DTRMV(U,N,N,   59, A,   64, X,     1)  MAXERR = .388799116737581D+11
Wrong RESULT  : DTRMV(U,N,N,   60, A,   64, X,     1)  MAXERR = .612787305842372D+11
Wrong RESULT  : DTRMV(U,N,N,   61, A,   64, X,     1)  MAXERR = .421813369205623D+11
Wrong RESULT  : DTRMV(U,N,N,   62, A,   64, X,     1)  MAXERR = .531064994597937D+11
Wrong RESULT  : DTRMV(U,N,N,   63, A,   64, X,     1)  MAXERR = .218143907859646D+12
Wrong RESULT  : DTRMV(U,N,N,   64, A,   64, X,     1)  MAXERR = .160600130355114D+12

The demonstration code is here: https://gist.github.com/grisuthedragon/32d8ad11e1d722414f921509b97d5507

[martin-frbg]

Is this with 0.2.20 or current git "develop" ?

[grisuthedragon]

With the current develop (b7cee00).

I checked the same on

• a SandyBride Xeon ( 16 Cores, E5-2690)
• a Westmere Xeon
• and my usual OpenPOWER 8 system

and the result ist the same for N = 1..16 the difference is in O(eps). For N=17..32 the difference is slightly increasing. For N>=32 is crashes completely as shown above. The fact on the OpenPower is interessting because the DTRMV routine is involved in some QR decomposition routines and so this may related to #1191

Changing the leading dimension does not influence the results as long N <= LDA <= 32.

If the number of threads is restricted to 1 everything is correct.

[martin-frbg]

Bisecting shows that this is not a recent failure, this behaviour appears to have been present in the last libGoto already.

[martin-frbg]

This code is indeed unchanged from (at least) GotoBLAS2-1.08 released in 2009.
The range calculations in trmv_thread.c and related level2 threading helpers have become a bit suspect in the context of #1089 already (where it appeared that the computed values could exceed the actual data size), but I do not really understand the workings of this code. (To be clear, undoing my PR #1262 does not help however)

[grisuthedragon]

Because of the fact that trmv is rarely used and from the computational density not the perfect one to have a threaded version I would turn of threading for it by setting nthreads=1 in interfaces/trmv.c. I tried this for the example above and the code where I got this error and it gives correct results without a noteworthy performance loss.

Furthermore, it seems to fix #1191 because in the example I used there, the dtpqrt function was used which has dtrmv inside its call graph. After setting the number of threads to one for the trmv routine the results stay correct for 100 independent runs.

[brada4]

Is it fast in many gigabytes scale too?

[grisuthedragon]

First thing must be to deliver the correct results and the the performance. Because it helps nobody to have incredible fast results, which are completely wrong.

[martin-frbg]

I would still prefer to understand just what is going wrong with your specific case, given that this bug appears to have gone unnoticed for years (and is not triggered by the compile-time tests).

[grisuthedragon]

The fact that the bug does not trigger the compile-time tests (as they are in NETLIB BLAS) is easy to explain. There are (together with the cblas interface) 735 calls to DTRMV and the largest parameter N is 9 and the largest leading dimension is 10 which is the area where my check code also gives the correct results. And the fact that it never get noticed. puuh.. From my experience I would say that it is rarely used and in most cases replaced by TRMM.

[martin-frbg]

Curioser and curioser... on my Kaby Lake laptop, repeated runs of your test code give varying patterns of "Correct" and "Wrong" results for all N beyond 16 - even with nthreads clamped to 1 in interface/trmv.c and confirming that only level2/trmv_U.c gets called. So there may be more than just one bug here. (Nehalem target behaves well with the nthreads=1 hack, but fails without it)

[martin-frbg]

Seems daxpy_microk_haswell-2.c is (subtly ?) broken on top of whatever is wrong with level2/trmv_thread.c - using either the generic C implementation or any of the older microkernels in its place brings the Haswell behaviour in line with Nehalem/Sandybridge/Atom

[martin-frbg]

Btw increasing N beyond 64 in your test code shows the results for even larger N still creeping away into madness, or is that simply a limitation in your code I am too dense to see now ?

[grisuthedragon]

I only did the example upto N=64 because that was the maximum size I need in my application. But in principal my code should work for larger N as well.

[martin-frbg]

Bad news then... Things apparently go wrong in trmv_U.c when (or shortly after) the problem size exceeds the DTB_ENTRIES (which is 64, the last "correct" calculation has N=67).
Update: a CORE2 build gets the right result (with nthreads forced to 1) up to much higher values of N,
probably by virtue of its DTB_ENTRIES count of 256)

[martin-frbg]

Unfortunately I am getting nowhere with this at the moment.