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
|
// 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"
"math"
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/lapack"
)
type Dgeconer interface {
Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
Dgetrier
Dlanger
}
func DgeconTest(t *testing.T, impl Dgeconer) {
rnd := rand.New(rand.NewSource(1))
for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
for _, lda := range []int{max(1, n), n + 3} {
dgeconTest(t, impl, rnd, n, lda)
}
}
}
func dgeconTest(t *testing.T, impl Dgeconer, rnd *rand.Rand, n, lda int) {
const ratioThresh = 10
// Generate a random square matrix A with elements uniformly in [-1,1).
a := make([]float64, max(0, (n-1)*lda+n))
for i := range a {
a[i] = 2*rnd.Float64() - 1
}
// Allocate work slices.
iwork := make([]int, n)
work := make([]float64, max(1, 4*n))
// Compute the LU factorization of A.
aFac := make([]float64, len(a))
copy(aFac, a)
ipiv := make([]int, n)
ok := impl.Dgetrf(n, n, aFac, lda, ipiv)
if !ok {
t.Fatalf("n=%v,lda=%v: bad matrix, Dgetrf failed", n, lda)
}
aFacCopy := make([]float64, len(aFac))
copy(aFacCopy, aFac)
// Compute the inverse A^{-1} from the LU factorization.
aInv := make([]float64, len(aFac))
copy(aInv, aFac)
ok = impl.Dgetri(n, aInv, lda, ipiv, work, len(work))
if !ok {
t.Fatalf("n=%v,lda=%v: bad matrix, Dgetri failed", n, lda)
}
for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
name := fmt.Sprintf("norm=%v,n=%v,lda=%v", string(norm), n, lda)
// Compute the norm of A and A^{-1}.
aNorm := impl.Dlange(norm, n, n, a, lda, work)
aInvNorm := impl.Dlange(norm, 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 the LU factorization and Dgecon.
rcondGot := impl.Dgecon(norm, n, aFac, lda, aNorm, work, iwork)
if !floats.Equal(aFac, aFacCopy) {
t.Errorf("%v: unexpected modification of aFac", 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)
}
// Check for corner-case values of anorm.
for _, anorm := range []float64{0, math.Inf(1), math.NaN()} {
rcondGot = impl.Dgecon(norm, n, aFac, lda, anorm, work, iwork)
if n == 0 {
if rcondGot != 1 {
t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=1", name, anorm, rcondGot)
}
continue
}
if math.IsNaN(anorm) {
if !math.IsNaN(rcondGot) {
t.Errorf("%v: NaN not propagated when anorm=NaN: got=%v", name, rcondGot)
}
continue
}
if rcondGot != 0 {
t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=0", name, anorm, rcondGot)
}
}
}
}
|
