1 // Copyright ©2017 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.
10 "golang.org/x/exp/rand"
12 "gonum.org/v1/gonum/blas"
13 "gonum.org/v1/gonum/blas/blas64"
14 "gonum.org/v1/gonum/floats"
15 "gonum.org/v1/gonum/lapack"
18 type Dtgsjaer interface {
20 Dtgsja(jobU, jobV, jobQ lapack.GSVDJob, m, p, n, k, l int, a []float64, lda int, b []float64, ldb int, tola, tolb float64, alpha, beta, u []float64, ldu int, v []float64, ldv int, q []float64, ldq int, work []float64) (cycles int, ok bool)
23 func DtgsjaTest(t *testing.T, impl Dtgsjaer) {
24 rnd := rand.New(rand.NewSource(1))
25 for cas, test := range []struct {
26 m, p, n, k, l, lda, ldb, ldu, ldv, ldq int
30 {m: 5, p: 5, n: 5, k: 2, l: 2, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
31 {m: 5, p: 5, n: 5, k: 4, l: 1, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
32 {m: 5, p: 5, n: 10, k: 2, l: 2, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
33 {m: 5, p: 5, n: 10, k: 4, l: 1, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
34 {m: 5, p: 5, n: 10, k: 4, l: 2, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
35 {m: 10, p: 5, n: 5, k: 2, l: 2, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
36 {m: 10, p: 5, n: 5, k: 4, l: 1, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
37 {m: 10, p: 10, n: 10, k: 5, l: 3, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
38 {m: 10, p: 10, n: 10, k: 6, l: 4, lda: 0, ldb: 0, ldu: 0, ldv: 0, ldq: 0, ok: true},
39 {m: 5, p: 5, n: 5, k: 2, l: 2, lda: 10, ldb: 10, ldu: 10, ldv: 10, ldq: 10, ok: true},
40 {m: 5, p: 5, n: 5, k: 4, l: 1, lda: 10, ldb: 10, ldu: 10, ldv: 10, ldq: 10, ok: true},
41 {m: 5, p: 5, n: 10, k: 2, l: 2, lda: 20, ldb: 20, ldu: 10, ldv: 10, ldq: 20, ok: true},
42 {m: 5, p: 5, n: 10, k: 4, l: 1, lda: 20, ldb: 20, ldu: 10, ldv: 10, ldq: 20, ok: true},
43 {m: 5, p: 5, n: 10, k: 4, l: 2, lda: 20, ldb: 20, ldu: 10, ldv: 10, ldq: 20, ok: true},
44 {m: 10, p: 5, n: 5, k: 2, l: 2, lda: 10, ldb: 10, ldu: 20, ldv: 10, ldq: 10, ok: true},
45 {m: 10, p: 5, n: 5, k: 4, l: 1, lda: 10, ldb: 10, ldu: 20, ldv: 10, ldq: 10, ok: true},
46 {m: 10, p: 10, n: 10, k: 5, l: 3, lda: 20, ldb: 20, ldu: 20, ldv: 20, ldq: 20, ok: true},
47 {m: 10, p: 10, n: 10, k: 6, l: 4, lda: 20, ldb: 20, ldu: 20, ldv: 20, ldq: 20, ok: true},
75 a := blockedUpperTriGeneral(m, n, k, l, lda, true, rnd)
76 aCopy := cloneGeneral(a)
77 b := blockedUpperTriGeneral(p, n, k, l, ldb, false, rnd)
78 bCopy := cloneGeneral(b)
80 tola := float64(max(m, n)) * impl.Dlange(lapack.NormFrob, m, n, a.Data, a.Stride, nil) * dlamchE
81 tolb := float64(max(p, n)) * impl.Dlange(lapack.NormFrob, p, n, b.Data, b.Stride, nil) * dlamchE
83 alpha := make([]float64, n)
84 beta := make([]float64, n)
86 work := make([]float64, 2*n)
88 u := nanGeneral(m, m, ldu)
89 v := nanGeneral(p, p, ldv)
90 q := nanGeneral(n, n, ldq)
92 _, ok := impl.Dtgsja(lapack.GSVDUnit, lapack.GSVDUnit, lapack.GSVDUnit,
105 t.Errorf("test %d unexpectedly did not converge", cas)
110 // Check orthogonality of U, V and Q.
111 if !isOrthonormal(u) {
112 t.Errorf("test %d: U is not orthogonal\n%+v", cas, u)
114 if !isOrthonormal(v) {
115 t.Errorf("test %d: V is not orthogonal\n%+v", cas, v)
117 if !isOrthonormal(q) {
118 t.Errorf("test %d: Q is not orthogonal\n%+v", cas, q)
121 // Check C^2 + S^2 = I.
122 var elements []float64
124 elements = alpha[k : k+l]
126 elements = alpha[k:m]
128 for i := range elements {
130 d := alpha[i]*alpha[i] + beta[i]*beta[i]
131 if !floats.EqualWithinAbsOrRel(d, 1, 1e-14, 1e-14) {
132 t.Errorf("test %d: alpha_%d^2 + beta_%d^2 != 1: got: %v", cas, i, i, d)
136 zeroR, d1, d2 := constructGSVDresults(n, p, m, k, l, a, b, alpha, beta)
138 // Check U^T*A*Q = D1*[ 0 R ].
139 uTmp := nanGeneral(m, n, n)
140 blas64.Gemm(blas.Trans, blas.NoTrans, 1, u, aCopy, 0, uTmp)
141 uAns := nanGeneral(m, n, n)
142 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, uTmp, q, 0, uAns)
144 d10r := nanGeneral(m, n, n)
145 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d1, zeroR, 0, d10r)
147 if !equalApproxGeneral(uAns, d10r, 1e-14) {
148 t.Errorf("test %d: U^T*A*Q != D1*[ 0 R ]\nU^T*A*Q:\n%+v\nD1*[ 0 R ]:\n%+v",
152 // Check V^T*B*Q = D2*[ 0 R ].
153 vTmp := nanGeneral(p, n, n)
154 blas64.Gemm(blas.Trans, blas.NoTrans, 1, v, bCopy, 0, vTmp)
155 vAns := nanGeneral(p, n, n)
156 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, vTmp, q, 0, vAns)
158 d20r := nanGeneral(p, n, n)
159 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d2, zeroR, 0, d20r)
161 if !equalApproxGeneral(vAns, d20r, 1e-14) {
162 t.Errorf("test %d: V^T*B*Q != D2*[ 0 R ]\nV^T*B*Q:\n%+v\nD2*[ 0 R ]:\n%+v",