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

   1 // Copyright ©2016 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         "fmt"
   9         "math"
  10         "testing"
  11
  12         "golang.org/x/exp/rand"
  13
  14         "gonum.org/v1/gonum/blas"
  15         "gonum.org/v1/gonum/blas/blas64"
  16 )
  17
  18 type Dorghrer interface {
  19         Dorghr(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int)
  20
  21         Dgehrder
  22 }
  23
  24 func DorghrTest(t *testing.T, impl Dorghrer) {
  25         rnd := rand.New(rand.NewSource(1))
  26
  27         for _, n := range []int{1, 2, 3, 4, 5, 6, 7, 8, 23, 34} {
  28                 for _, extra := range []int{0, 1, 13} {
  29                         for _, optwork := range []bool{true, false} {
  30                                 for cas := 0; cas < 100; cas++ {
  31                                         ilo := rnd.Intn(n)
  32                                         ihi := rnd.Intn(n)
  33                                         if ilo > ihi {
  34                                                 ilo, ihi = ihi, ilo
  35                                         }
  36                                         testDorghr(t, impl, n, ilo, ihi, extra, optwork, rnd)
  37                                 }
  38                         }
  39                 }
  40         }
  41         testDorghr(t, impl, 0, 0, -1, 0, false, rnd)
  42         testDorghr(t, impl, 0, 0, -1, 0, true, rnd)
  43 }
  44
  45 func testDorghr(t *testing.T, impl Dorghrer, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
  46         const tol = 1e-14
  47
  48         // Construct the matrix A with elementary reflectors and scalar factors tau.
  49         a := randomGeneral(n, n, n+extra, rnd)
  50         var tau []float64
  51         if n > 1 {
  52                 tau = nanSlice(n - 1)
  53         }
  54         work := nanSlice(max(1, n)) // Minimum work for Dgehrd.
  55         impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
  56
  57         // Extract Q for later comparison.
  58         q := eye(n, n)
  59         qCopy := cloneGeneral(q)
  60         for j := ilo; j < ihi; j++ {
  61                 h := eye(n, n)
  62                 v := blas64.Vector{
  63                         Inc:  1,
  64                         Data: make([]float64, n),
  65                 }
  66                 v.Data[j+1] = 1
  67                 for i := j + 2; i < ihi+1; i++ {
  68                         v.Data[i] = a.Data[i*a.Stride+j]
  69                 }
  70                 blas64.Ger(-tau[j], v, v, h)
  71                 copy(qCopy.Data, q.Data)
  72                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
  73         }
  74
  75         if optwork {
  76                 work = nanSlice(1)
  77                 impl.Dorghr(n, ilo, ihi, a.Data, a.Stride, tau, work, -1)
  78                 work = nanSlice(int(work[0]))
  79         } else {
  80                 work = nanSlice(max(1, ihi-ilo))
  81         }
  82         impl.Dorghr(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
  83
  84         prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v, optwork=%v", n, ilo, ihi, extra, optwork)
  85         if !generalOutsideAllNaN(a) {
  86                 t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
  87         }
  88         if !isOrthonormal(a) {
  89                 t.Errorf("%v: A is not orthogonal\n%v", prefix, a.Data)
  90         }
  91         for i := 0; i < n; i++ {
  92                 for j := 0; j < n; j++ {
  93                         aij := a.Data[i*a.Stride+j]
  94                         qij := q.Data[i*q.Stride+j]
  95                         if math.Abs(aij-qij) > tol {
  96                                 t.Errorf("%v: unexpected value of A[%v,%v]. want %v, got %v", prefix, i, j, qij, aij)
  97                         }
  98                 }
  99         }
 100 }