diff --git a/lapack/testlapack/dtrexc.go b/lapack/testlapack/dtrexc.go index 4b20675e..59d0d3d0 100644 --- a/lapack/testlapack/dtrexc.go +++ b/lapack/testlapack/dtrexc.go @@ -27,7 +27,7 @@ func DtrexcTest(t *testing.T, impl Dtrexcer) { for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 18, 31, 53} { for _, extra := range []int{0, 1, 11} { for cas := 0; cas < 100; cas++ { - tmat := randomSchurCanonical(n, n+extra, rnd) + tmat, _, _ := randomSchurCanonical(n, n+extra, rnd) ifst := rnd.Intn(n) ilst := rnd.Intn(n) testDtrexc(t, impl, compq, tmat, ifst, ilst, extra, rnd) @@ -38,7 +38,7 @@ func DtrexcTest(t *testing.T, impl Dtrexcer) { for _, compq := range []lapack.UpdateSchurComp{lapack.UpdateSchurNone, lapack.UpdateSchur} { for _, extra := range []int{0, 1, 11} { - tmat := randomSchurCanonical(0, extra, rnd) + tmat, _, _ := randomSchurCanonical(0, extra, rnd) testDtrexc(t, impl, compq, tmat, 0, 0, extra, rnd) } } diff --git a/lapack/testlapack/general.go b/lapack/testlapack/general.go index 075a4006..96a16f54 100644 --- a/lapack/testlapack/general.go +++ b/lapack/testlapack/general.go @@ -131,8 +131,8 @@ func randomHessenberg(n, stride int, rnd *rand.Rand) blas64.General { // form, that is, block upper triangular with 1×1 and 2×2 diagonal blocks where // each 2×2 diagonal block has its diagonal elements equal and its off-diagonal // elements of opposite sign. -func randomSchurCanonical(n, stride int, rnd *rand.Rand) blas64.General { - t := randomGeneral(n, n, stride, rnd) +func randomSchurCanonical(n, stride int, rnd *rand.Rand) (t blas64.General, wr, wi []float64) { + t = randomGeneral(n, n, stride, rnd) // Zero out the lower triangle. for i := 0; i < t.Rows; i++ { for j := 0; j < i; j++ { @@ -141,23 +141,49 @@ func randomSchurCanonical(n, stride int, rnd *rand.Rand) blas64.General { } // Randomly create 2×2 diagonal blocks. for i := 0; i < t.Rows; { - if i == t.Rows-1 || rnd.Float64() < 0.5 { + p := rnd.Float64() + var r float64 + switch { + case p < 0.5: + // A half of real parts of eigenvalues will be "normal". + r = t.Data[i*t.Stride+i] + case p < 0.75: + // A quarter will be very small. + r = dlamchS + default: + // A quarter will be zero. + } + // Store the actual value back at the diagonal of T. + t.Data[i*t.Stride+i] = r + + // A half of eigenvalues will be real. + if rnd.Float64() < 0.5 || i == t.Rows-1 { // 1×1 block. + wr = append(wr, r) + wi = append(wi, 0) i++ continue } + // 2×2 block. // Diagonal elements equal. - t.Data[(i+1)*t.Stride+i+1] = t.Data[i*t.Stride+i] + t.Data[(i+1)*t.Stride+i+1] = r + // Element under diagonal very small or "normal". + c := dlamchS + if rnd.Float64() < 0.5 { + c = rnd.NormFloat64() + } // Off-diagonal elements of opposite sign. - c := rnd.NormFloat64() if math.Signbit(c) == math.Signbit(t.Data[i*t.Stride+i+1]) { c *= -1 } t.Data[(i+1)*t.Stride+i] = c + wr = append(wr, r, r) + c = math.Sqrt(math.Abs(t.Data[i*t.Stride+i+1])) * math.Sqrt(math.Abs(c)) + wi = append(wi, c, -c) i += 2 } - return t + return t, wr, wi } // blockedUpperTriGeneral returns a normal random, general matrix in the form