vendor/gonum.org/v1/gonum/lapack/testlapack/dgehd2.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 Dgehd2er interface {
  19         Dgehd2(n, ilo, ihi int, a []float64, lda int, tau, work []float64)
  20 }
  21
  22 func Dgehd2Test(t *testing.T, impl Dgehd2er) {
  23         rnd := rand.New(rand.NewSource(1))
  24         for _, n := range []int{1, 2, 3, 4, 5, 7, 10, 30} {
  25                 for _, extra := range []int{0, 1, 13} {
  26                         for cas := 0; cas < 100; cas++ {
  27                                 testDgehd2(t, impl, n, extra, rnd)
  28                         }
  29                 }
  30         }
  31 }
  32
  33 func testDgehd2(t *testing.T, impl Dgehd2er, n, extra int, rnd *rand.Rand) {
  34         ilo := rnd.Intn(n)
  35         ihi := rnd.Intn(n)
  36         if ilo > ihi {
  37                 ilo, ihi = ihi, ilo
  38         }
  39
  40         tau := nanSlice(n - 1)
  41         work := nanSlice(n)
  42
  43         a := randomGeneral(n, n, n+extra, rnd)
  44         // NaN out elements under the diagonal except
  45         // for the [ilo:ihi,ilo:ihi] block.
  46         for i := 1; i <= ihi; i++ {
  47                 for j := 0; j < min(ilo, i); j++ {
  48                         a.Data[i*a.Stride+j] = math.NaN()
  49                 }
  50         }
  51         for i := ihi + 1; i < n; i++ {
  52                 for j := 0; j < i; j++ {
  53                         a.Data[i*a.Stride+j] = math.NaN()
  54                 }
  55         }
  56         aCopy := a
  57         aCopy.Data = make([]float64, len(a.Data))
  58         copy(aCopy.Data, a.Data)
  59
  60         impl.Dgehd2(n, ilo, ihi, a.Data, a.Stride, tau, work)
  61
  62         prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v", n, ilo, ihi, extra)
  63
  64         // Check any invalid modifications of a.
  65         if !generalOutsideAllNaN(a) {
  66                 t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
  67         }
  68         for i := ilo; i <= ihi; i++ {
  69                 for j := 0; j < min(ilo, i); j++ {
  70                         if !math.IsNaN(a.Data[i*a.Stride+j]) {
  71                                 t.Errorf("%v: expected NaN at A[%v,%v]", prefix, i, j)
  72                         }
  73                 }
  74         }
  75         for i := ihi + 1; i < n; i++ {
  76                 for j := 0; j < i; j++ {
  77                         if !math.IsNaN(a.Data[i*a.Stride+j]) {
  78                                 t.Errorf("%v: expected NaN at A[%v,%v]", prefix, i, j)
  79                         }
  80                 }
  81         }
  82         for i := 0; i <= ilo; i++ {
  83                 for j := i; j < ilo+1; j++ {
  84                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
  85                                 t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
  86                         }
  87                 }
  88                 for j := ihi + 1; j < n; j++ {
  89                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
  90                                 t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
  91                         }
  92                 }
  93         }
  94         for i := ihi + 1; i < n; i++ {
  95                 for j := i; j < n; j++ {
  96                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
  97                                 t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
  98                         }
  99                 }
 100         }
 101
 102         // Check that tau has been assigned properly.
 103         for i, v := range tau {
 104                 if i < ilo || i >= ihi {
 105                         if !math.IsNaN(v) {
 106                                 t.Errorf("%v: expected NaN at tau[%v]", prefix, i)
 107                         }
 108                 } else {
 109                         if math.IsNaN(v) {
 110                                 t.Errorf("%v: unexpected NaN at tau[%v]", prefix, i)
 111                         }
 112                 }
 113         }
 114
 115         // Extract Q and check that it is orthogonal.
 116         q := blas64.General{
 117                 Rows:   n,
 118                 Cols:   n,
 119                 Stride: n,
 120                 Data:   make([]float64, n*n),
 121         }
 122         for i := 0; i < q.Rows; i++ {
 123                 q.Data[i*q.Stride+i] = 1
 124         }
 125         qCopy := q
 126         qCopy.Data = make([]float64, len(q.Data))
 127         for j := ilo; j < ihi; j++ {
 128                 h := blas64.General{
 129                         Rows:   n,
 130                         Cols:   n,
 131                         Stride: n,
 132                         Data:   make([]float64, n*n),
 133                 }
 134                 for i := 0; i < h.Rows; i++ {
 135                         h.Data[i*h.Stride+i] = 1
 136                 }
 137                 v := blas64.Vector{
 138                         Inc:  1,
 139                         Data: make([]float64, n),
 140                 }
 141                 v.Data[j+1] = 1
 142                 for i := j + 2; i < ihi+1; i++ {
 143                         v.Data[i] = a.Data[i*a.Stride+j]
 144                 }
 145                 blas64.Ger(-tau[j], v, v, h)
 146                 copy(qCopy.Data, q.Data)
 147                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
 148         }
 149         if !isOrthonormal(q) {
 150                 t.Errorf("%v: Q is not orthogonal\nQ=%v", prefix, q)
 151         }
 152
 153         // Overwrite NaN elements of aCopy with zeros
 154         // (we will multiply with it below).
 155         for i := 1; i <= ihi; i++ {
 156                 for j := 0; j < min(ilo, i); j++ {
 157                         aCopy.Data[i*aCopy.Stride+j] = 0
 158                 }
 159         }
 160         for i := ihi + 1; i < n; i++ {
 161                 for j := 0; j < i; j++ {
 162                         aCopy.Data[i*aCopy.Stride+j] = 0
 163                 }
 164         }
 165
 166         // Construct Q^T * AOrig * Q and check that it is
 167         // equal to A from Dgehd2.
 168         aq := blas64.General{
 169                 Rows:   n,
 170                 Cols:   n,
 171                 Stride: n,
 172                 Data:   make([]float64, n*n),
 173         }
 174         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aCopy, q, 0, aq)
 175         qaq := blas64.General{
 176                 Rows:   n,
 177                 Cols:   n,
 178                 Stride: n,
 179                 Data:   make([]float64, n*n),
 180         }
 181         blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, aq, 0, qaq)
 182         for i := ilo; i <= ihi; i++ {
 183                 for j := ilo; j <= ihi; j++ {
 184                         qaqij := qaq.Data[i*qaq.Stride+j]
 185                         if j < i-1 {
 186                                 if math.Abs(qaqij) > 1e-14 {
 187                                         t.Errorf("%v: Q^T*A*Q is not upper Hessenberg, [%v,%v]=%v", prefix, i, j, qaqij)
 188                                 }
 189                                 continue
 190                         }
 191                         diff := qaqij - a.Data[i*a.Stride+j]
 192                         if math.Abs(diff) > 1e-14 {
 193                                 t.Errorf("%v: Q^T*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)
 194                         }
 195                 }
 196         }
 197 }