diff --git a/testlapack/dgebal.go b/testlapack/dgebal.go index 52beedc7..911ccb72 100644 --- a/testlapack/dgebal.go +++ b/testlapack/dgebal.go @@ -6,7 +6,6 @@ package testlapack import ( "fmt" - "math" "math/rand" "testing" @@ -24,15 +23,9 @@ func DgebalTest(t *testing.T, impl Dgebaler) { for _, job := range []lapack.Job{lapack.None, lapack.Permute, lapack.Scale, lapack.PermuteScale} { for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 10, 18, 31, 53, 100} { - var pNonzero float64 - if n == 0 || n == 1 { - pNonzero = 0.5 - } else { - pNonzero = math.Ceil(0.02*float64(n)) / float64(n) - } for _, extra := range []int{0, 11} { for cas := 0; cas < 100; cas++ { - a := randomSparseGeneral(n, n, n+extra, pNonzero, rnd) + a := unbalancedSparseGeneral(n, n, n+extra, 2*n, rnd) testDgebal(t, impl, job, a) } } @@ -57,6 +50,10 @@ func testDgebal(t *testing.T, impl Dgebaler, job lapack.Job, a blas64.General) { prefix := fmt.Sprintf("Case job=%v, n=%v, extra=%v", job, n, extra) + if !generalOutsideAllNaN(a) { + t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data) + } + if n == 0 { if ilo != 0 { t.Errorf("%v: unexpected ilo when n=0. Want 0, got %v", prefix, n, ilo) diff --git a/testlapack/general.go b/testlapack/general.go index a740d541..4d4458bd 100644 --- a/testlapack/general.go +++ b/testlapack/general.go @@ -1014,17 +1014,18 @@ func isUpperTriangular(a blas64.General) bool { return true } -// randomSparseGeneral returns an m×n general matrix whose each element is set -// to a non-zero value with probability p. -func randomSparseGeneral(m, n, stride int, p float64, rnd *rand.Rand) blas64.General { - a := nanGeneral(m, n, stride) - for i := 0; i < m; i++ { - for j := 0; j < n; j++ { - if rnd.NormFloat64() < p { - a.Data[i*stride+j] = rnd.NormFloat64() - } else { - a.Data[i*stride+j] = 0 - } +// unbalancedSparseGeneral returns an m×n dense matrix with a random sparse +// structure consisting of nz nonzero elements. The matrix will be unbalanced by +// multiplying each element randomly by its row or column index. +func unbalancedSparseGeneral(m, n, stride int, nonzeros int, rnd *rand.Rand) blas64.General { + a := zeros(m, n, stride) + for k := 0; k < nonzeros; k++ { + i := rnd.Intn(n) + j := rnd.Intn(n) + if rnd.Float64() < 0.5 { + a.Data[i*stride+j] = float64(i+1) * rnd.NormFloat64() + } else { + a.Data[i*stride+j] = float64(j+1) * rnd.NormFloat64() } } return a