some question about rotmg

[wxcstc]

Well,when i used rotmg function,the same input had different result between cublas and openblas.
input is as follows:
d1=5.9e-8
d2=5.960464e-8
x1=1.0
y1=150.0
param[5]={0.0, 0.0, 0.0, 0.0, 0.0}

And the output of openblas:
d1=0.999956
d2=0.989812
x1=0.036623
y1=150.0
param[5]={-1.0, 0.000002, -0.000244, 0.000244, 0.000002}

And the output of cublas:
d1=0.999956
d2=0.989812
x1=0.036623
y1=150.0
param[5]={-1.0, 0.000002, -0.000244, 1.0, 0.000002}

I referred some papers in order to know the implementation of the function.And I found that in some situation, this function will have a scaling operation.It seems that cublas does not have this operation? It is a bug of cublas?I think it's a bug of cublas, but i'm not sure.

[martin-frbg]

Ahem, did you notice the reply to your identical message in the nvidia cuda forum where the moderator asked you for a standalone test case to replicate the issue ?

[martin-frbg]

On the other hand, a testcase with your parameters linked against the netlib reference BLAS does give the same result as you got with cublas, so you may be onto something.

[martin-frbg]

The changed result comes from 025fc91 in May 2014 which must have been #365. There, the problem was repeated execution of a code segment based on the flag value set by its first invocation, now the added conditional to prevent this also seems to block a valid action in the first pass. Perhaps what is needed is a second flag, or perhaps I fail to understand the logic - this function was rewritten from old fortran-style code that jumped back and forth, neither the old nor the new version are easy to read.
Pinging @btracey as he suggested the fix that seems to cause problems now.

[martin-frbg]

I further notice that ATLAS - which was given as reference in #365 - seems to avoid setting any flags within while() loops, and indeed the scaling there is done in dedicated inner loops as needed, rather than looping over the entire logic...

[brada4]

Sorry, but are you using SROTMG or DROTMG? With former the inputs look marginal towards precision and may overflow something. There is trivial geometry behind this function. can you try netlib BLAS Like one built by your linux distro for ultimate truth?

[martin-frbg]

@brada4 drotmg, I already compared with netlib above (which agrees with cublas but not openblas). The supposed "fix" from 365 that broke it looks more like a bad workaround for a fundamental problem: IMHO rotmg should have "if" where it has "while" and do the while loop for scaling inside that "if" after reacting to (and resetting) dflag. What it does now is potentially flip dflag around with every iteration of the rescaling, which simply does not make sense.
I will do a PR tomorrow or over the weekend as time permits.

[brada4]

Yes, it looks like manual f2c from reference function...

[wxcstc]

@martin-frbg @brada4 I don't think it's the reason of dflag in while loop. In my understanding, if have multiple cycles，in the first cycle，dflag and other variables will be changed. And in next cycle，dflag has been -1，it will not go into the branch which to set dflag and other variables. So，using if or while to set dflag and other variables dosen't affect the final result. How do you think?

[martin-frbg]

The code for in next cycle, dflag has been -1, it will not go into the branch was added in #365 to prevent this unwanted side effect of cycling. The problem with it is that in your case, dflag is -1 from the start and the added code completely prevents it from changing the variables accordingly. My change separates the scaling cycle from the dflag setting code.

[wxcstc]

ok，I have already understood how this result occurred. Now I have another confusion，why should judge dflag two times？Can I extract the two judgments and put them in front of if(*dd1 != ZERO)？I am not very clear about the principle.With reference to this paper "basic linear algebra subprograms for fortran usage" http://www.cs.utexas.edu/users/kincaid/blas.pdf, I can understand most of the code.But still some do not understand. I did not find the statement to support the operation of judging the flag. Maybe my english is poor 😔. Can you answer my doubts？Thank you very much~😀

[kortschak]

I don't think this is completely fixed.

In attempting to fix our implementation I find that the whole situation is SNAFU; to get my head around it, I have done a literal translation of the FORTRAN77 code here rendered as this Go code. This fails to agree with my translation of the current NETLIB implementation here rendered as this Go code under our test suite:

--- FAIL: TestDrotmg (0.00s)
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Big_RD2_Big_Flag_0: expected -1 [4096 -3584 1792 4096], found -1 [4096 -4096 1 4096]
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Big_RD2_Big_Flag_1: expected -1 [2340.5714285714284 -4096 4096 4681.142857142857], found -1 [2340.5714285714284 -4096 1 4681.142857142857]
	level1double.go:2007: drotmg/drotmg_f77 H mismatch D1_Big_D2_Small_Flag_1: expected -1 [2.8671999999999997e-26 -0.000244140625 4096 2.44140625e-16], found -1 [2.8671999999999997e-26 -0.000244140625 1 2.44140625e-16]
	level1double.go:2007: drotmg/drotmg_f77 H mismatch D1_Small_D2_Big_Flag_1: expected -1 [2.3731773997569866e+10 -1.6777216e+07 0.000244140625 1.6777216e-07], found -1 [2.3731773997569866e+10 -4096 1 1.6777216e-07]
	level1double.go:2007: drotmg/drotmg_f77 H mismatch OpenBLAS#1452: expected -1 [1.6110934624105326e-06 -0.000244140625 0.000244140625 1.6276041666666668e-06], found -1 [1.6110934624105326e-06 -0.000244140625 1 1.6276041666666668e-06]

This makes the whole situation fairly confusing, given that neither of these actually satisfy the documentation's claims AFAICS (though the documentation is unclear and seem to be incorrect - vis the reference to DY2 which does not exists - though using DY1 does not help).

The situation is worse with the OpenBLAS implementation used as a CGO back-end for Go tests (some of these appear to be simple scaling differences, but others not):

--- FAIL: TestDrotmg (0.00s)
	level1double.go:1994: drotmg/drotmg_f77 rd1 mismatch RD1_Big_RD2_Big_Flag_0: expected 68.96627824858757, found 1.1570621468926554e+09
	level1double.go:1997: drotmg/drotmg_f77 rd2 mismatch RD1_Big_RD2_Big_Flag_0: expected 34.483139124293785, found 5.785310734463277e+08
	level1double.go:2000: drotmg/drotmg_f77 rx1 mismatch RD1_Big_RD2_Big_Flag_0: expected 45312, found 11.0625
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Big_RD2_Big_Flag_0: expected -1 [4096 -3584 1792 4096], found -1 [1 -1 1 1]
	level1double.go:1994: drotmg/drotmg_f77 rd1 mismatch RD1_Big_RD2_Big_Flag_1: expected 57.6914092640818, found 9.679012345679014e+08
	level1double.go:1997: drotmg/drotmg_f77 rd2 mismatch RD1_Big_RD2_Big_Flag_1: expected 28.8457046320409, found 4.839506172839507e+08
	level1double.go:2000: drotmg/drotmg_f77 rx1 mismatch RD1_Big_RD2_Big_Flag_1: expected 47396.57142857142, found 11.57142857142857
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Big_RD2_Big_Flag_1: expected -1 [2340.5714285714284 -4096 4096 4681.142857142857], found -1 [0.5714285714285714 -1 1 1.1428571428571428]
	level1double.go:1994: drotmg/drotmg_f77 rd1 mismatch RD1_Big_RD2_Med_Flag_0: expected 1.1920927762985347, found 1.9999998000000197e+07
	level1double.go:2000: drotmg/drotmg_f77 rx1 mismatch RD1_Big_RD2_Med_Flag_0: expected 32768.0032768, found 8.0000008
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Big_RD2_Med_Flag_0: expected -1 [4096 -1 0.0004096 1], found -1 [1 -1 1e-07 1]
	level1double.go:1994: drotmg/drotmg_f77 rd1 mismatch RD1_Big_RD2_Med_Flag_1: expected 1192.0928955078125, found 2e+10
	level1double.go:2000: drotmg/drotmg_f77 rx1 mismatch RD1_Big_RD2_Med_Flag_1: expected 3.2768e+14, found 8e+10
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Big_RD2_Med_Flag_1: expected -1 [4.096e-17 -1 4096 1e-10], found -1 [1e-20 -1 1 1e-10]
	level1double.go:1994: drotmg/drotmg_f77 rd1 mismatch D1_Big_D2_Small_Flag_1: expected 119.20928955078125, found 2e+09
	level1double.go:2000: drotmg/drotmg_f77 rx1 mismatch D1_Big_D2_Small_Flag_1: expected 3.2768e+10, found 8e+06
	level1double.go:2007: drotmg/drotmg_f77 H mismatch D1_Big_D2_Small_Flag_1: expected -1 [2.8671999999999997e-26 -0.000244140625 4096 2.44140625e-16], found -1 [6.999999999999999e-30 -0.000244140625 1 2.44140625e-16]
	level1double.go:1997: drotmg/drotmg_f77 rd2 mismatch RD1_Med_RD2_Big_Flag_0: expected 1191.9736981379988, found 1.9998000199980003e+10
	level1double.go:2007: drotmg/drotmg_f77 H mismatch RD1_Med_RD2_Big_Flag_0: expected -1 [1 -0.0004096 1000 4096], found -1 [1 -1e-07 1000 1]
	level1double.go:1997: drotmg/drotmg_f77 rd2 mismatch D1_Med_D2_Big_Flag_1: expected 1192.092835903171, found 1.9999999000000053e+10
	level1double.go:2007: drotmg/drotmg_f77 H mismatch D1_Med_D2_Big_Flag_1: expected -1 [50 -4096 1 4.096e-06], found -1 [50 -1 1 1e-09]
	level1double.go:1997: drotmg/drotmg_f77 rd2 mismatch D1_Small_D2_Big_Flag_1: expected 216.1836123957717, found 6.085027726432532e+16
	level1double.go:2007: drotmg/drotmg_f77 H mismatch D1_Small_D2_Big_Flag_1: expected -1 [2.3731773997569866e+10 -1.6777216e+07 0.000244140625 1.6777216e-07], found -1 [2.3731773997569866e+10 -1 1 1e-14]
	level1double.go:2007: drotmg/drotmg_f77 H mismatch OpenBLAS#1452: expected -1 [1.6110934624105326e-06 -0.000244140625 0.000244140625 1.6276041666666668e-06], found -1 [1.6110934624105326e-06 -0.000244140625 1 1.6276041666666668e-06]

Note that in all cases the mismatched H is in the DFLAG=-1 state, so all elements must satisfy the check.