1 // Copyright ©2015 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
11 "golang.org/x/exp/rand"
13 "gonum.org/v1/gonum/blas"
14 "gonum.org/v1/gonum/floats"
15 "gonum.org/v1/gonum/lapack"
18 type Dtrconer interface {
20 Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64
23 func DtrconTest(t *testing.T, impl Dtrconer) {
24 rnd := rand.New(rand.NewSource(1))
25 // Hand crafted tests.
26 for _, test := range []struct {
55 condOne: 0.137704918032787,
56 condInf: 0.157894736842105,
79 condOne: 0.157894736842105,
80 condInf: 0.137704918032787,
84 work := make([]float64, 3*test.n)
86 work[i] = rnd.Float64()
88 iwork := make([]int, test.n)
89 for i := range iwork {
90 iwork[i] = int(rnd.Int31())
92 aCopy := make([]float64, len(test.a))
94 condOne := impl.Dtrcon(lapack.MaxColumnSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork)
95 if math.Abs(condOne-test.condOne) > 1e-14 {
96 t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne)
98 if !floats.Equal(aCopy, test.a) {
99 t.Errorf("a modified during call")
101 condInf := impl.Dtrcon(lapack.MaxRowSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork)
102 if math.Abs(condInf-test.condInf) > 1e-14 {
103 t.Errorf("Inf norm mismatch. Want %v, got %v.", test.condInf, condInf)
105 if !floats.Equal(aCopy, test.a) {
106 t.Errorf("a modified during call")
110 // Dtrcon does not match the Dgecon output in many cases. See
111 // https://github.com/xianyi/OpenBLAS/issues/636
112 // TODO(btracey): Uncomment this when the mismatch between Dgecon and Dtrcon
115 // Randomized tests against Dgecon.
116 for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
117 for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
118 for _, test := range []struct {
124 for trial := 0; trial < 1; trial++ {
130 a := make([]float64, n*lda)
133 a[i] = float64(i + 2)
137 a[i] = rnd.NormFloat64()
141 aDense := make([]float64, len(a))
142 if uplo == blas.Upper {
143 for i := 0; i < n; i++ {
144 for j := i; j < n; j++ {
145 aDense[i*lda+j] = a[i*lda+j]
149 for i := 0; i < n; i++ {
150 for j := 0; j <= i; j++ {
151 aDense[i*lda+j] = a[i*lda+j]
155 if diag == blas.Unit {
156 for i := 0; i < n; i++ {
161 ipiv := make([]int, n)
162 work := make([]float64, 4*n)
163 denseOne := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work)
164 denseInf := impl.Dlange(lapack.MaxRowSum, n, n, aDense, lda, work)
166 aDenseLU := make([]float64, len(aDense))
167 copy(aDenseLU, aDense)
168 impl.Dgetrf(n, n, aDenseLU, lda, ipiv)
169 iwork := make([]int, n)
170 want := impl.Dgecon(lapack.MaxColumnSum, n, aDenseLU, lda, denseOne, work, iwork)
171 got := impl.Dtrcon(lapack.MaxColumnSum, uplo, diag, n, a, lda, work, iwork)
172 if math.Abs(want-got) > 1e-14 {
173 t.Errorf("One norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got)
175 want = impl.Dgecon(lapack.MaxRowSum, n, aDenseLU, lda, denseInf, work, iwork)
176 got = impl.Dtrcon(lapack.MaxRowSum, uplo, diag, n, a, lda, work, iwork)
177 if math.Abs(want-got) > 1e-14 {
178 t.Errorf("Inf norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got)