lapack/testlapack: rework test for Dtrcon
MODIFY
lapack/testlapack/dtrcon.go
lapack/testlapack/dtrcon.go
package testlapack
import (
- "math"
+ "fmt"
"testing"
"golang.org/x/exp/rand"
@@ -16,170 +16,116 @@ import (
)
type Dtrconer interface {
- Dgeconer
Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64
+
+ Dtrtri(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int) bool
+ Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a []float64, lda int, work []float64) float64
}
func DtrconTest(t *testing.T, impl Dtrconer) {
rnd := rand.New(rand.NewSource(1))
- // Hand crafted tests.
- for _, test := range []struct {
- a []float64
- n int
- uplo blas.Uplo
- diag blas.Diag
- condOne float64
- condInf float64
- }{
- {
- a: []float64{
- 8, 5, 6,
- 0, 7, 8,
- 0, 0, 6,
- },
- n: 3,
- uplo: blas.Upper,
- diag: blas.Unit,
- condOne: 1.0 / 645,
- condInf: 1.0 / 480,
- },
- {
- a: []float64{
- 8, 5, 6,
- 0, 7, 8,
- 0, 0, 6,
- },
- n: 3,
- uplo: blas.Upper,
- diag: blas.NonUnit,
- condOne: 0.137704918032787,
- condInf: 0.157894736842105,
- },
- {
- a: []float64{
- 8, 0, 0,
- 5, 7, 0,
- 6, 8, 6,
- },
- n: 3,
- uplo: blas.Lower,
- diag: blas.Unit,
- condOne: 1.0 / 480,
- condInf: 1.0 / 645,
- },
- {
- a: []float64{
- 8, 0, 0,
- 5, 7, 0,
- 6, 8, 6,
- },
- n: 3,
- uplo: blas.Lower,
- diag: blas.NonUnit,
- condOne: 0.157894736842105,
- condInf: 0.137704918032787,
- },
- } {
- lda := test.n
- work := make([]float64, 3*test.n)
- for i := range work {
- work[i] = rnd.Float64()
- }
- iwork := make([]int, test.n)
- for i := range iwork {
- iwork[i] = int(rnd.Int31())
- }
- aCopy := make([]float64, len(test.a))
- copy(aCopy, test.a)
- condOne := impl.Dtrcon(lapack.MaxColumnSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork)
- if math.Abs(condOne-test.condOne) > 1e-14 {
- t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne)
- }
- if !floats.Equal(aCopy, test.a) {
- t.Errorf("a modified during call")
- }
- condInf := impl.Dtrcon(lapack.MaxRowSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork)
- if math.Abs(condInf-test.condInf) > 1e-14 {
- t.Errorf("Inf norm mismatch. Want %v, got %v.", test.condInf, condInf)
- }
- if !floats.Equal(aCopy, test.a) {
- t.Errorf("a modified during call")
- }
- }
-
- // Dtrcon does not match the Dgecon output in many cases. See
- // https://github.com/xianyi/OpenBLAS/issues/636
- // TODO(btracey): Uncomment this when the mismatch between Dgecon and Dtrcon
- // is understood.
- /*
- // Randomized tests against Dgecon.
+ for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
- for _, test := range []struct {
- n, lda int
- }{
- {3, 0},
- {4, 9},
- } {
- for trial := 0; trial < 1; trial++ {
- n := test.n
- lda := test.lda
- if lda == 0 {
- lda = n
- }
- a := make([]float64, n*lda)
- if trial == 0 {
- for i := range a {
- a[i] = float64(i + 2)
+ for _, lda := range []int{max(1, n), n + 3} {
+ for _, mattype := range []int{0, 1, 2} {
+ dtrconTest(t, impl, rnd, uplo, diag, n, lda, mattype)
+ }
+ }
+ }
}
- } else {
- for i := range a {
- a[i] = rnd.NormFloat64()
}
}
- aDense := make([]float64, len(a))
+func dtrconTest(t *testing.T, impl Dtrconer, rnd *rand.Rand, uplo blas.Uplo, diag blas.Diag, n, lda, mattype int) {
+ const ratioThresh = 10
+
+ a := make([]float64, max(0, (n-1)*lda+n))
+ for i := range a {
+ a[i] = rnd.Float64()
+ }
+ switch mattype {
+ default:
+ panic("bad mattype")
+ case 0:
+ // Matrix filled with consecutive integer values.
+ // For lapack.MaxRowSum norm (infinity-norm) these matrices
+ // sometimes lead to a slightly inaccurate estimate of the condition
+ // number.
+ c := 2.0
+ for i := 0; i < n; i++ {
+ for j := 0; j < n; j++ {
+ a[i*lda+j] = c
+ c += 1
+ }
+ }
+ case 1:
+ // Identity matrix.
if uplo == blas.Upper {
for i := 0; i < n; i++ {
- for j := i; j < n; j++ {
- aDense[i*lda+j] = a[i*lda+j]
+ for j := i + 1; j < n; j++ {
+ a[i*lda+j] = 0
}
}
} else {
for i := 0; i < n; i++ {
- for j := 0; j <= i; j++ {
- aDense[i*lda+j] = a[i*lda+j]
+ for j := 0; j < i; j++ {
+ a[i*lda+j] = 0
}
}
}
- if diag == blas.Unit {
+ if diag == blas.NonUnit {
for i := 0; i < n; i++ {
- aDense[i*lda+i] = 1
+ a[i*lda+i] = 1
}
}
+ case 2:
+ // Matrix filled with random values uniformly in [-1,1).
+ // These matrices often lead to a slightly inaccurate estimate
+ // of the condition number.
+ for i := 0; i < n; i++ {
+ for j := 0; j < n; j++ {
+ a[i*lda+j] = 2*rnd.Float64() - 1
+ }
+ }
+ }
+ aCopy := make([]float64, len(a))
+ copy(aCopy, a)
- ipiv := make([]int, n)
- work := make([]float64, 4*n)
- denseOne := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work)
- denseInf := impl.Dlange(lapack.MaxRowSum, n, n, aDense, lda, work)
+ // Compute the inverse A^{-1}.
+ aInv := make([]float64, len(a))
+ copy(aInv, a)
+ ok := impl.Dtrtri(uplo, diag, n, aInv, lda)
+ if !ok {
+ t.Fatalf("uplo=%v,diag=%v,n=%v,lda=%v,mattype=%v: bad matrix, Dtrtri failed", string(uplo), string(diag), n, lda, mattype)
+ }
- aDenseLU := make([]float64, len(aDense))
- copy(aDenseLU, aDense)
- impl.Dgetrf(n, n, aDenseLU, lda, ipiv)
+ work := make([]float64, 3*n)
iwork := make([]int, n)
- want := impl.Dgecon(lapack.MaxColumnSum, n, aDenseLU, lda, denseOne, work, iwork)
- got := impl.Dtrcon(lapack.MaxColumnSum, uplo, diag, n, a, lda, work, iwork)
- if math.Abs(want-got) > 1e-14 {
- t.Errorf("One norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got)
- }
- want = impl.Dgecon(lapack.MaxRowSum, n, aDenseLU, lda, denseInf, work, iwork)
- got = impl.Dtrcon(lapack.MaxRowSum, uplo, diag, n, a, lda, work, iwork)
- if math.Abs(want-got) > 1e-14 {
- t.Errorf("Inf norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got)
- }
+ for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
+ name := fmt.Sprintf("norm=%v,uplo=%v,diag=%v,n=%v,lda=%v,mattype=%v", string(norm), string(uplo), string(diag), n, lda, mattype)
+
+ // Compute the norm of A and A^{-1}.
+ aNorm := impl.Dlantr(norm, uplo, diag, n, n, a, lda, work)
+ aInvNorm := impl.Dlantr(norm, uplo, diag, n, n, aInv, lda, work)
+
+ // Compute a good estimate of the condition number
+ // rcondWant := 1/(norm(A) * norm(inv(A)))
+ rcondWant := 1.0
+ if aNorm > 0 && aInvNorm > 0 {
+ rcondWant = 1 / aNorm / aInvNorm
}
+
+ // Compute an estimate of rcond using Dtrcon.
+ rcondGot := impl.Dtrcon(norm, uplo, diag, n, a, lda, work, iwork)
+ if !floats.Equal(a, aCopy) {
+ t.Errorf("%v: unexpected modification of a", name)
}
+
+ ratio := rCondTestRatio(rcondGot, rcondWant)
+ if ratio >= ratioThresh {
+ t.Errorf("%v: unexpected value of rcond; got=%v, want=%v (ratio=%v)",
+ name, rcondGot, rcondWant, ratio)
}
}
- */
}

MODIFY
lapack/testlapack/dtrtri.go
lapack/testlapack/dtrtri.go
@@ -14,7 +14,6 @@ import (
)
type Dtrtrier interface {
- Dtrconer
Dtrtri(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int) bool
}

lapack/testlapack: rework test for Dtrcon