diff --git a/testlapack/dlaexc.go b/testlapack/dlaexc.go index 3e695f5..240b6b4 100644 --- a/testlapack/dlaexc.go +++ b/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,6 +79,7 @@ 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) } @@ -72,6 +87,9 @@ func testDlaexc(t *testing.T, impl Dlaexcer, wantq bool, n, j1, n1, n2, extra in 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) } }