vendor/gonum.org/v1/gonum/lapack/testlapack/dormhr.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 Dormhrer interface {
  19         Dormhr(side blas.Side, trans blas.Transpose, m, n, ilo, ihi int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
  20
  21         Dgehrder
  22 }
  23
  24 func DormhrTest(t *testing.T, impl Dormhrer) {
  25         rnd := rand.New(rand.NewSource(1))
  26
  27         for _, side := range []blas.Side{blas.Left, blas.Right} {
  28                 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
  29                         for _, m := range []int{1, 2, 3, 4, 5, 8, 9, 10, 23} {
  30                                 for _, n := range []int{1, 2, 3, 4, 5, 8, 9, 10, 23} {
  31                                         for _, extra := range []int{0, 1, 13} {
  32                                                 for cas := 0; cas < 10; cas++ {
  33                                                         nq := m
  34                                                         if side == blas.Right {
  35                                                                 nq = n
  36                                                         }
  37                                                         ilo := rnd.Intn(nq)
  38                                                         ihi := rnd.Intn(nq)
  39                                                         if ilo > ihi {
  40                                                                 ilo, ihi = ihi, ilo
  41                                                         }
  42                                                         testDormhr(t, impl, side, trans, m, n, ilo, ihi, extra, true, rnd)
  43                                                         testDormhr(t, impl, side, trans, m, n, ilo, ihi, extra, false, rnd)
  44                                                 }
  45                                         }
  46                                 }
  47                         }
  48                 }
  49         }
  50         for _, side := range []blas.Side{blas.Left, blas.Right} {
  51                 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
  52                         testDormhr(t, impl, side, trans, 0, 0, 0, -1, 0, true, rnd)
  53                         testDormhr(t, impl, side, trans, 0, 0, 0, -1, 0, false, rnd)
  54                 }
  55         }
  56 }
  57
  58 func testDormhr(t *testing.T, impl Dormhrer, side blas.Side, trans blas.Transpose, m, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
  59         const tol = 1e-14
  60
  61         var nq, nw int
  62         switch side {
  63         case blas.Left:
  64                 nq = m
  65                 nw = n
  66         case blas.Right:
  67                 nq = n
  68                 nw = m
  69         }
  70
  71         // Compute the elementary reflectors and tau.
  72         a := randomGeneral(nq, nq, nq+extra, rnd)
  73         var tau []float64
  74         if nq > 1 {
  75                 tau = nanSlice(nq - 1)
  76         }
  77         work := nanSlice(max(1, nq)) // Minimum work for Dgehrd.
  78         impl.Dgehrd(nq, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
  79
  80         // Construct Q from the elementary reflectors in a and from tau.
  81         q := eye(nq, nq)
  82         qCopy := eye(nq, nq)
  83         for j := ilo; j < ihi; j++ {
  84                 h := eye(nq, nq)
  85                 v := blas64.Vector{
  86                         Inc:  1,
  87                         Data: make([]float64, nq),
  88                 }
  89                 v.Data[j+1] = 1
  90                 for i := j + 2; i < ihi+1; i++ {
  91                         v.Data[i] = a.Data[i*a.Stride+j]
  92                 }
  93                 blas64.Ger(-tau[j], v, v, h)
  94                 copy(qCopy.Data, q.Data)
  95                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
  96         }
  97
  98         c := randomGeneral(m, n, n+extra, rnd)
  99
 100         // Compute the product of Q and C explicitly.
 101         qc := randomGeneral(m, n, n+extra, rnd)
 102         if side == blas.Left {
 103                 blas64.Gemm(trans, blas.NoTrans, 1, q, c, 0, qc)
 104         } else {
 105                 blas64.Gemm(blas.NoTrans, trans, 1, c, q, 0, qc)
 106         }
 107
 108         // Compute the product of Q and C using Dormhr.
 109         if optwork {
 110                 work = nanSlice(1)
 111                 impl.Dormhr(side, trans, m, n, ilo, ihi, nil, a.Stride, nil, nil, c.Stride, work, -1)
 112                 work = nanSlice(int(work[0]))
 113         } else {
 114                 work = nanSlice(max(1, nw))
 115         }
 116         impl.Dormhr(side, trans, m, n, ilo, ihi, a.Data, a.Stride, tau, c.Data, c.Stride, work, len(work))
 117
 118         // Compare the two answers.
 119         prefix := fmt.Sprintf("Case side=%v, trans=%v, m=%v, n=%v, ilo=%v, ihi=%v, extra=%v, optwork=%v",
 120                 side, trans, m, n, ilo, ihi, extra, optwork)
 121         if !generalOutsideAllNaN(c) {
 122                 t.Errorf("%v: out-of-range write to C\n%v", prefix, c.Data)
 123         }
 124         for i := 0; i < m; i++ {
 125                 for j := 0; j < n; j++ {
 126                         cij := c.Data[i*c.Stride+j]
 127                         qcij := qc.Data[i*qc.Stride+j]
 128                         if math.Abs(cij-qcij) > tol {
 129                                 t.Errorf("%v: unexpected value of the QC product at [%v,%v]: want %v, got %v", prefix, i, j, qcij, cij)
 130                         }
 131                 }
 132         }
 133 }