testlapack: fix and update test for Dlaexc

testlapack/dlaexc.go

@@ -16,6 +16,8 @@ import (
 
 type Dlaexcer interface {
 	Dlaexc(wantq bool, n int, t []float64, ldt int, q []float64, ldq int, j1, n1, n2 int, work []float64) bool
+
+	Dlanv2(a, b, c, d float64) (aa, bb, cc, dd float64, rt1r, rt1i, rt2r, rt2i float64, cs, sn float64)
 }
 
 func DlaexcTest(t *testing.T, impl Dlaexcer) {
@@ -47,12 +49,24 @@ func testDlaexc(t *testing.T, impl Dlaexcer, wantq bool, n, j1, n1, n2, extra in
 	}
 	// Make any 2x2 diagonal block to be in Schur canonical form.
 	if n1 == 2 {
+		// Diagonal elements equal.
 		tmat.Data[(j1+1)*tmat.Stride+j1+1] = tmat.Data[j1*tmat.Stride+j1]
-		tmat.Data[(j1+1)*tmat.Stride+j1] = tmat.Data[j1*tmat.Stride+j1+1]
+		// Off-diagonal elements of opposite sign.
+		c := rnd.NormFloat64()
+		if math.Signbit(c) == math.Signbit(tmat.Data[j1*tmat.Stride+j1+1]) {
+			c *= -1
+		}
+		tmat.Data[(j1+1)*tmat.Stride+j1] = c
 	}
 	if n2 == 2 {
+		// Diagonal elements equal.
 		tmat.Data[(j1+n1+1)*tmat.Stride+j1+n1+1] = tmat.Data[(j1+n1)*tmat.Stride+j1+n1]
-		tmat.Data[(j1+n1+1)*tmat.Stride+j1+n1] = tmat.Data[(j1+n1)*tmat.Stride+j1+n1+1]
+		// Off-diagonal elements of opposite sign.
+		c := rnd.NormFloat64()
+		if math.Signbit(c) == math.Signbit(tmat.Data[(j1+n1)*tmat.Stride+j1+n1+1]) {
+			c *= -1
+		}
+		tmat.Data[(j1+n1+1)*tmat.Stride+j1+n1] = c
 	}
 	tmatCopy := cloneGeneral(tmat)
 	var q, qCopy blas64.General
@@ -65,13 +79,17 @@ func testDlaexc(t *testing.T, impl Dlaexcer, wantq bool, n, j1, n1, n2, extra in
 	ok := impl.Dlaexc(wantq, n, tmat.Data, tmat.Stride, q.Data, q.Stride, j1, n1, n2, work)
 
 	prefix := fmt.Sprintf("Case n=%v, j1=%v, n1=%v, n2=%v, wantq=%v, extra=%v", n, j1, n1, n2, wantq, extra)
+
 	if !generalOutsideAllNaN(tmat) {
 		t.Errorf("%v: out-of-range write to T\n", prefix)
 	}
 	if wantq && !generalOutsideAllNaN(q) {
 		t.Errorf("%v: out-of-range write to Q\n", prefix)
 	}
 
+	if !ok {
+		t.Logf("%v: Dlaexc returned false")
+	}
 	if !ok || n1 == 0 || n2 == 0 || j1+n1 >= n {
 		// Check that T is not modified.
 		for i := 0; i < n; i++ {
@@ -110,23 +128,74 @@ func testDlaexc(t *testing.T, impl Dlaexcer, wantq bool, n, j1, n1, n2, extra in
 			}
 		}
 	}
-	// Check that T is modified at the diagonal as expected.
-	for i := 0; i < n1; i++ {
-		for j := 0; j < n1; j++ {
-			got := tmat.Data[(j1+i)*tmat.Stride+j1+j]
-			want := tmatCopy.Data[(j1+n1+i)*tmatCopy.Stride+j1+n1+j]
-			if want != got {
-				t.Errorf("%v: unexpected value of T[%v,%v]. want %v, got %v\n", prefix, j1+i, j1+j, want, got)
-			}
+	if n1 == 1 {
+		// 1×1 blocks are swapped exactly.
+		got := tmat.Data[(j1+n2)*tmat.Stride+j1+n2]
+		want := tmatCopy.Data[j1*tmatCopy.Stride+j1]
+		if want != got {
+			t.Errorf("%v: unexpected value of T[%v,%v]. want %v, got %v\n", prefix, j1+n2, j1+n2, want, got)
+		}
+	} else {
+		// Check that the swapped 2×2 block is in canonical Schur form.
+		// The n1×n1 block is now located at T[j1+n2,j1+n2].
+		a := tmat.Data[(j1+n2)*tmat.Stride+j1+n2]
+		b := tmat.Data[(j1+n2)*tmat.Stride+j1+n2+1]
+		c := tmat.Data[(j1+n2+1)*tmat.Stride+j1+n2]
+		d := tmat.Data[(j1+n2+1)*tmat.Stride+j1+n2+1]
+		if !isCanonicalSchur(a, b, c, d) {
+			t.Errorf("%v: 2×2 block at T[%v,%v] not in canonical Schur form", prefix, j1+n2, j1+n2)
+		}
+
+		// Check that the swapped 2×2 block has the same eigenvalues.
+		_, _, _, _, rt1rGot, rt1iGot, rt2rGot, rt2iGot, _, _ := impl.Dlanv2(a, b, c, d)
+		// The n1×n1 block was originally located at T[j1,j1].
+		a = tmatCopy.Data[j1*tmatCopy.Stride+j1]
+		b = tmatCopy.Data[j1*tmatCopy.Stride+j1+1]
+		c = tmatCopy.Data[(j1+1)*tmatCopy.Stride+j1]
+		d = tmatCopy.Data[(j1+1)*tmatCopy.Stride+j1+1]
+		_, _, _, _, rt1rWant, rt1iWant, rt2rWant, rt2iWant, _, _ := impl.Dlanv2(a, b, c, d)
+		if math.Abs(rt1rWant-rt1rGot) > tol || math.Abs(rt1iWant-rt1iGot) > tol {
+			t.Errorf("%v: unexpected first eigenvalue of 2×2 block at T[%v,%v]. Want (%v,%v), got (%v,%v)",
+				prefix, j1+n2, j1+n2, rt1rWant, rt1iWant, rt1rGot, rt1iGot)
+		}
+		if math.Abs(rt2rWant-rt2rGot) > tol || math.Abs(rt2iWant-rt2iGot) > tol {
+			t.Errorf("%v: unexpected second eigenvalue of 2×2 block at T[%v,%v]. Want (%v,%v), got (%v,%v)",
+				prefix, j1+n2, j1+n2, rt2rWant, rt2iWant, rt2rGot, rt2iGot)
 		}
 	}
-	for i := 0; i < n2; i++ {
-		for j := 0; j < n2; j++ {
-			got := tmat.Data[(j1+n1+i)*tmat.Stride+j1+n1+j]
-			want := tmatCopy.Data[(j1+i)*tmatCopy.Stride+j1+j]
-			if want != got {
-				t.Errorf("%v: unexpected value of T[%v,%v]. want %v, got %v\n", prefix, j1+n1+i, j1+n1+j, want, got)
-			}
+	if n2 == 1 {
+		// 1×1 blocks are swapped exactly.
+		got := tmat.Data[j1*tmat.Stride+j1]
+		want := tmatCopy.Data[(j1+n1)*tmatCopy.Stride+j1+n1]
+		if want != got {
+			t.Errorf("%v: unexpected value of T[%v,%v]. want %v, got %v\n", prefix, j1, j1, want, got)
+		}
+	} else {
+		// Check that the swapped 2×2 block is in canonical Schur form.
+		// The n2×n2 block is now located at T[j1,j1].
+		a := tmat.Data[j1*tmat.Stride+j1]
+		b := tmat.Data[j1*tmat.Stride+j1+1]
+		c := tmat.Data[(j1+1)*tmat.Stride+j1]
+		d := tmat.Data[(j1+1)*tmat.Stride+j1+1]
+		if !isCanonicalSchur(a, b, c, d) {
+			t.Errorf("%v: 2×2 block at T[%v,%v] not in canonical Schur form", prefix, j1, j1)
+		}
+
+		// Check that the swapped 2×2 block has the same eigenvalues.
+		_, _, _, _, rt1rGot, rt1iGot, rt2rGot, rt2iGot, _, _ := impl.Dlanv2(a, b, c, d)
+		// The n2×n2 block was originally located at T[j1+n1,j1+n1].
+		a = tmatCopy.Data[(j1+n1)*tmatCopy.Stride+j1+n1]
+		b = tmatCopy.Data[(j1+n1)*tmatCopy.Stride+j1+n1+1]
+		c = tmatCopy.Data[(j1+n1+1)*tmatCopy.Stride+j1+n1]
+		d = tmatCopy.Data[(j1+n1+1)*tmatCopy.Stride+j1+n1+1]
+		_, _, _, _, rt1rWant, rt1iWant, rt2rWant, rt2iWant, _, _ := impl.Dlanv2(a, b, c, d)
+		if math.Abs(rt1rWant-rt1rGot) > tol || math.Abs(rt1iWant-rt1iGot) > tol {
+			t.Errorf("%v: unexpected first eigenvalue of 2×2 block at T[%v,%v]. Want (%v,%v), got (%v,%v)",
+				prefix, j1, j1, rt1rWant, rt1iWant, rt1rGot, rt1iGot)
+		}
+		if math.Abs(rt2rWant-rt2rGot) > tol || math.Abs(rt2iWant-rt2iGot) > tol {
+			t.Errorf("%v: unexpected second eigenvalue of 2×2 block at T[%v,%v]. Want (%v,%v), got (%v,%v)",
+				prefix, j1, j1, rt2rWant, rt2iWant, rt2rGot, rt2iGot)
 		}
 	}