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
|
// Copyright ©2017 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 (
"math"
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
"gonum.org/v1/gonum/floats/scalar"
)
type Dlags2er interface {
Dlags2(upper bool, a1, a2, a3, b1, b2, b3 float64) (csu, snu, csv, snv, csq, snq float64)
}
func Dlags2Test(t *testing.T, impl Dlags2er) {
rnd := rand.New(rand.NewSource(1))
for _, upper := range []bool{true, false} {
for i := 0; i < 100; i++ {
// Generate randomly the elements of a 2×2 matrix A
// [ a1 a2 ] or [ a1 0 ]
// [ 0 a3 ] [ a2 a3 ]
a1 := rnd.Float64()
a2 := rnd.Float64()
a3 := rnd.Float64()
// Generate randomly the elements of a 2×2 matrix B.
// [ b1 b2 ] or [ b1 0 ]
// [ 0 b3 ] [ b2 b3 ]
b1 := rnd.Float64()
b2 := rnd.Float64()
b3 := rnd.Float64()
// Compute orthogonal matrices U, V, Q
// U = [ csu snu ], V = [ csv snv ], Q = [ csq snq ]
// [ -snu csu ] [ -snv csv ] [ -snq csq ]
// that transform A and B.
csu, snu, csv, snv, csq, snq := impl.Dlags2(upper, a1, a2, a3, b1, b2, b3)
// Check that U, V, Q are orthogonal matrices (their
// determinant is equal to 1).
detU := det2x2(csu, snu, -snu, csu)
if !scalar.EqualWithinAbsOrRel(math.Abs(detU), 1, 1e-14, 1e-14) {
t.Errorf("U not orthogonal: det(U)=%v", detU)
}
detV := det2x2(csv, snv, -snv, csv)
if !scalar.EqualWithinAbsOrRel(math.Abs(detV), 1, 1e-14, 1e-14) {
t.Errorf("V not orthogonal: det(V)=%v", detV)
}
detQ := det2x2(csq, snq, -snq, csq)
if !scalar.EqualWithinAbsOrRel(math.Abs(detQ), 1, 1e-14, 1e-14) {
t.Errorf("Q not orthogonal: det(Q)=%v", detQ)
}
// Create U, V, Q explicitly as dense matrices.
u := blas64.General{
Rows: 2,
Cols: 2,
Stride: 2,
Data: []float64{csu, snu, -snu, csu},
}
v := blas64.General{
Rows: 2,
Cols: 2,
Stride: 2,
Data: []float64{csv, snv, -snv, csv},
}
q := blas64.General{
Rows: 2,
Cols: 2,
Stride: 2,
Data: []float64{csq, snq, -snq, csq},
}
// Create A and B explicitly as dense matrices.
a := blas64.General{Rows: 2, Cols: 2, Stride: 2}
b := blas64.General{Rows: 2, Cols: 2, Stride: 2}
if upper {
a.Data = []float64{a1, a2, 0, a3}
b.Data = []float64{b1, b2, 0, b3}
} else {
a.Data = []float64{a1, 0, a2, a3}
b.Data = []float64{b1, 0, b2, b3}
}
tmp := blas64.General{Rows: 2, Cols: 2, Stride: 2, Data: make([]float64, 4)}
// Transform A as Uᵀ*A*Q.
blas64.Gemm(blas.Trans, blas.NoTrans, 1, u, a, 0, tmp)
blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, q, 0, a)
// Transform B as Vᵀ*A*Q.
blas64.Gemm(blas.Trans, blas.NoTrans, 1, v, b, 0, tmp)
blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, q, 0, b)
// Extract elements of transformed A and B that should be equal to zero.
var gotA, gotB float64
if upper {
gotA = a.Data[1]
gotB = b.Data[1]
} else {
gotA = a.Data[2]
gotB = b.Data[2]
}
// Check that they are indeed zero.
if !scalar.EqualWithinAbsOrRel(gotA, 0, 1e-14, 1e-14) {
t.Errorf("unexpected non-zero value for zero triangle of Uᵀ*A*Q: %v", gotA)
}
if !scalar.EqualWithinAbsOrRel(gotB, 0, 1e-14, 1e-14) {
t.Errorf("unexpected non-zero value for zero triangle of Vᵀ*B*Q: %v", gotB)
}
}
}
}
func det2x2(a, b, c, d float64) float64 { return a*d - b*c }
|
