vendor/gonum.org/v1/gonum/lapack/testlapack/dtgsja.go

   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.
   4
   5 package testlapack
   6
   7 import (
   8         "testing"
   9
  10         "golang.org/x/exp/rand"
  11
  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"
  16 )
  17
  18 type Dtgsjaer interface {
  19         Dlanger
  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)
  21 }
  22
  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
  27
  28                 ok bool
  29         }{
  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},
  48         } {
  49                 m := test.m
  50                 p := test.p
  51                 n := test.n
  52                 k := test.k
  53                 l := test.l
  54                 lda := test.lda
  55                 if lda == 0 {
  56                         lda = n
  57                 }
  58                 ldb := test.ldb
  59                 if ldb == 0 {
  60                         ldb = n
  61                 }
  62                 ldu := test.ldu
  63                 if ldu == 0 {
  64                         ldu = m
  65                 }
  66                 ldv := test.ldv
  67                 if ldv == 0 {
  68                         ldv = p
  69                 }
  70                 ldq := test.ldq
  71                 if ldq == 0 {
  72                         ldq = n
  73                 }
  74
  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)
  79
  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
  82
  83                 alpha := make([]float64, n)
  84                 beta := make([]float64, n)
  85
  86                 work := make([]float64, 2*n)
  87
  88                 u := nanGeneral(m, m, ldu)
  89                 v := nanGeneral(p, p, ldv)
  90                 q := nanGeneral(n, n, ldq)
  91
  92                 _, ok := impl.Dtgsja(lapack.GSVDUnit, lapack.GSVDUnit, lapack.GSVDUnit,
  93                         m, p, n, k, l,
  94                         a.Data, a.Stride,
  95                         b.Data, b.Stride,
  96                         tola, tolb,
  97                         alpha, beta,
  98                         u.Data, u.Stride,
  99                         v.Data, v.Stride,
 100                         q.Data, q.Stride,
 101                         work)
 102
 103                 if !ok {
 104                         if test.ok {
 105                                 t.Errorf("test %d unexpectedly did not converge", cas)
 106                         }
 107                         continue
 108                 }
 109
 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)
 113                 }
 114                 if !isOrthonormal(v) {
 115                         t.Errorf("test %d: V is not orthogonal\n%+v", cas, v)
 116                 }
 117                 if !isOrthonormal(q) {
 118                         t.Errorf("test %d: Q is not orthogonal\n%+v", cas, q)
 119                 }
 120
 121                 // Check C^2 + S^2 = I.
 122                 var elements []float64
 123                 if m-k-l >= 0 {
 124                         elements = alpha[k : k+l]
 125                 } else {
 126                         elements = alpha[k:m]
 127                 }
 128                 for i := range elements {
 129                         i += k
 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)
 133                         }
 134                 }
 135
 136                 zeroR, d1, d2 := constructGSVDresults(n, p, m, k, l, a, b, alpha, beta)
 137
 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)
 143
 144                 d10r := nanGeneral(m, n, n)
 145                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d1, zeroR, 0, d10r)
 146
 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",
 149                                 cas, uAns, d10r)
 150                 }
 151
 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)
 157
 158                 d20r := nanGeneral(p, n, n)
 159                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d2, zeroR, 0, d20r)
 160
 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",
 163                                 cas, vAns, d20r)
 164                 }
 165         }
 166 }