diff --git a/lapack/testlapack/dlartg.go b/lapack/testlapack/dlartg.go index e0c414a1..c2c82cfa 100644 --- a/lapack/testlapack/dlartg.go +++ b/lapack/testlapack/dlartg.go @@ -20,24 +20,32 @@ type Dlartger interface { func DlartgTest(t *testing.T, impl Dlartger) { const tol = 1e-14 // safmn2 and safmx2 are copied from native.Dlartg. + // safmn2 ~ 2*10^{-146} safmn2 := math.Pow(dlamchB, math.Trunc(math.Log(dlamchS/dlamchE)/math.Log(dlamchB)/2)) + // safmx2 ~ 5*10^145 safmx2 := 1 / safmn2 rnd := rand.New(rand.NewSource(1)) for i := 0; i < 1000; i++ { + // Generate randomly huge, tiny, and "normal" input arguments to Dlartg. var f float64 var fHuge bool switch rnd.Intn(3) { case 0: // Huge f. fHuge = true - f = math.Pow(10, 10-20*rnd.Float64()) * safmx2 + // scale is in the range (10^{-10}, 10^10]. + scale := math.Pow(10, 10-20*rnd.Float64()) + // f is in the range (5*10^135, 5*10^155]. + f = scale * safmx2 case 1: // Tiny f. + // f is in the range (2*10^{-156}, 2*10^{-136}]. f = math.Pow(10, 10-20*rnd.Float64()) * safmn2 default: f = rnd.NormFloat64() } if rnd.Intn(2) == 0 { + // Sometimes change the sign of f. f *= -1 } @@ -58,26 +66,38 @@ func DlartgTest(t *testing.T, impl Dlartger) { g *= -1 } + // Generate a plane rotation so that + // [ cs sn] * [f] = [r] + // [-sn cs] [g] = [0] cs, sn, r := impl.Dlartg(f, g) + // Check that the first equation holds. rWant := cs*f + sn*g if !floats.EqualWithinAbsOrRel(math.Abs(rWant), math.Abs(r), tol, tol) { t.Errorf("Case f=%v,g=%v: unexpected r. Want %v, got %v", f, g, rWant, r) } + // Check that cs and sn define a plane rotation. The 2×2 matrix + // has orthogonal columns by construction, so only check that + // the columns/rows have unit norm. oneTest := cs*cs + sn*sn if math.Abs(oneTest-1) > tol { t.Errorf("Case f=%v,g=%v: expected cs^2+sn^2==1, got %v", f, g, oneTest) } if !fHuge && !gHuge { + // Check that the second equation holds. + // If both numbers are huge, cancellation errors make + // this test unreliable. zeroTest := -sn*f + cs*g if math.Abs(zeroTest) > tol { t.Errorf("Case f=%v,g=%v: expected zero, got %v", f, g, zeroTest) } } + // Check that cs is positive as documented. if math.Abs(f) > math.Abs(g) && cs < 0 { t.Errorf("Case f=%v,g=%v: unexpected negative cs %v", f, g, cs) } } + // Check other documented special cases. for i := 0; i < 100; i++ { cs, sn, _ := impl.Dlartg(rnd.NormFloat64(), 0) if cs != 1 {