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

   1 // Copyright ©2015 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         "math"
   9         "testing"
  10
  11         "golang.org/x/exp/rand"
  12
  13         "gonum.org/v1/gonum/blas"
  14         "gonum.org/v1/gonum/blas/blas64"
  15 )
  16
  17 type Dlasq1er interface {
  18         Dlasq1(n int, d, e, work []float64) int
  19         Dgetrfer
  20 }
  21
  22 func Dlasq1Test(t *testing.T, impl Dlasq1er) {
  23         rnd := rand.New(rand.NewSource(1))
  24         bi := blas64.Implementation()
  25         // TODO(btracey): Increase the size of this test when we have a more numerically
  26         // stable way to test the singular values.
  27         for _, n := range []int{1, 2, 5, 8} {
  28                 work := make([]float64, 4*n)
  29                 d := make([]float64, n)
  30                 e := make([]float64, n-1)
  31                 for cas := 0; cas < 1; cas++ {
  32                         for i := range work {
  33                                 work[i] = rnd.Float64()
  34                         }
  35                         for i := range d {
  36                                 d[i] = rnd.NormFloat64() + 10
  37                         }
  38                         for i := range e {
  39                                 e[i] = rnd.NormFloat64()
  40                         }
  41                         ldm := n
  42                         m := make([]float64, n*ldm)
  43                         // Set up the matrix
  44                         for i := 0; i < n; i++ {
  45                                 m[i*ldm+i] = d[i]
  46                                 if i != n-1 {
  47                                         m[(i+1)*ldm+i] = e[i]
  48                                 }
  49                         }
  50
  51                         ldmm := n
  52                         mm := make([]float64, n*ldmm)
  53                         bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, m, ldm, m, ldm, 0, mm, ldmm)
  54
  55                         impl.Dlasq1(n, d, e, work)
  56
  57                         // Check that they are singular values. The
  58                         // singular values are the square roots of the
  59                         // eigenvalues of X^T * X
  60                         mmCopy := make([]float64, len(mm))
  61                         copy(mmCopy, mm)
  62                         ipiv := make([]int, n)
  63                         for elem, sv := range d[0:n] {
  64                                 copy(mm, mmCopy)
  65                                 lambda := sv * sv
  66                                 for i := 0; i < n; i++ {
  67                                         mm[i*ldm+i] -= lambda
  68                                 }
  69
  70                                 // Compute LU.
  71                                 ok := impl.Dgetrf(n, n, mm, ldmm, ipiv)
  72                                 if !ok {
  73                                         // Definitely singular.
  74                                         continue
  75                                 }
  76                                 // Compute determinant
  77                                 var logdet float64
  78                                 for i := 0; i < n; i++ {
  79                                         v := mm[i*ldm+i]
  80                                         logdet += math.Log(math.Abs(v))
  81                                 }
  82                                 if math.Exp(logdet) > 2 {
  83                                         t.Errorf("Incorrect singular value. n = %d, cas = %d, elem = %d, det = %v", n, cas, elem, math.Exp(logdet))
  84                                 }
  85                         }
  86                 }
  87         }
  88 }