1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
// Copyright ©2015 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package testlapack
import (
"fmt"
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/lapack"
)
type Dtrconer interface {
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))
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 _, lda := range []int{max(1, n), n + 3} {
for _, mattype := range []int{0, 1, 2} {
dtrconTest(t, impl, rnd, uplo, diag, n, lda, mattype)
}
}
}
}
}
}
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 + 1; j < n; j++ {
a[i*lda+j] = 0
}
}
} else {
for i := 0; i < n; i++ {
for j := 0; j < i; j++ {
a[i*lda+j] = 0
}
}
}
if diag == blas.NonUnit {
for i := 0; i < n; i++ {
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)
// 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)
}
work := make([]float64, 3*n)
iwork := make([]int, n)
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)
}
}
}
|
