vendor/gonum.org/v1/gonum/lapack/internal/testdata/dlaqr5test/main.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 // This program generates test data for Dlaqr5. Test cases are stored in
   6 // gzip-compressed JSON file testlapack/testdata/dlaqr5data.json.gz which is
   7 // read during testing by testlapack/dlaqr5.go.
   8 //
   9 // This program uses cgo to call Fortran version of DLAQR5. Therefore, matrices
  10 // passed to the Fortran routine are in column-major format but are written into
  11 // the output file in row-major format.
  12 package main
  13
  14 import (
  15         "compress/gzip"
  16         "encoding/json"
  17         "log"
  18         "os"
  19         "path/filepath"
  20
  21         "golang.org/x/exp/rand"
  22
  23         "gonum.org/v1/gonum/lapack/internal/testdata/netlib"
  24 )
  25
  26 type Dlaqr5Test struct {
  27         WantT          bool
  28         N              int
  29         NShifts        int
  30         KTop, KBot     int
  31         ShiftR, ShiftI []float64
  32         H              []float64
  33
  34         HWant []float64
  35         ZWant []float64
  36 }
  37
  38 func main() {
  39         file, err := os.Create(filepath.FromSlash("../../../testlapack/testdata/dlaqr5data.json.gz"))
  40         if err != nil {
  41                 log.Fatal(err)
  42         }
  43         defer file.Close()
  44         w := gzip.NewWriter(file)
  45
  46         rnd := rand.New(rand.NewSource(1))
  47
  48         var tests []Dlaqr5Test
  49         for _, wantt := range []bool{true, false} {
  50                 for _, n := range []int{2, 3, 4, 5, 6, 7, 11} {
  51                         for k := 0; k <= min(5, n); k++ {
  52                                 npairs := k
  53                                 if npairs == 0 {
  54                                         npairs = 2 * n
  55                                 }
  56                                 for ktop := 0; ktop < n-1; ktop++ {
  57                                         for kbot := ktop + 1; kbot < n; kbot++ {
  58                                                 sr, si := shiftpairs(npairs, rnd)
  59                                                 nshfts := len(sr)
  60
  61                                                 v := genrand(nshfts/2, 3, rnd)
  62                                                 u := genrand(3*nshfts-3, 3*nshfts-3, rnd)
  63                                                 wh := genrand(3*nshfts-3, n, rnd)
  64                                                 nh := n
  65                                                 wv := genrand(n, 3*nshfts-3, rnd)
  66                                                 nv := n
  67
  68                                                 h := hessrand(n, rnd)
  69                                                 if ktop > 0 {
  70                                                         h[ktop+(ktop-1)*n] = 0
  71                                                 }
  72                                                 if kbot < n-1 {
  73                                                         h[kbot+1+kbot*n] = 0
  74                                                 }
  75                                                 hin := make([]float64, len(h))
  76                                                 copy(hin, h)
  77                                                 z := eye(n)
  78
  79                                                 netlib.Dlaqr5(wantt, true, 2,
  80                                                         n, ktop+1, kbot+1,
  81                                                         nshfts, sr, si,
  82                                                         h, n,
  83                                                         1, n, z, n,
  84                                                         v, 3,
  85                                                         u, 3*nshfts-3,
  86                                                         nh, wh, nh,
  87                                                         nv, wv, 3*nshfts-3)
  88
  89                                                 tests = append(tests, Dlaqr5Test{
  90                                                         WantT:   wantt,
  91                                                         N:       n,
  92                                                         NShifts: nshfts,
  93                                                         KTop:    ktop,
  94                                                         KBot:    kbot,
  95                                                         ShiftR:  sr,
  96                                                         ShiftI:  si,
  97                                                         H:       rowMajor(n, n, hin),
  98                                                         HWant:   rowMajor(n, n, h),
  99                                                         ZWant:   rowMajor(n, n, z),
 100                                                 })
 101                                         }
 102                                 }
 103                         }
 104                 }
 105         }
 106         json.NewEncoder(w).Encode(tests)
 107
 108         err = w.Close()
 109         if err != nil {
 110                 log.Fatal(err)
 111         }
 112 }
 113
 114 // genrand returns a general r×c matrix with random entries.
 115 func genrand(r, c int, rnd *rand.Rand) []float64 {
 116         m := make([]float64, r*c)
 117         for i := range m {
 118                 m[i] = rnd.NormFloat64()
 119         }
 120         return m
 121 }
 122
 123 // eye returns an identity matrix of order n.
 124 func eye(n int) []float64 {
 125         m := make([]float64, n*n)
 126         for i := 0; i < n*n; i += n + 1 {
 127                 m[i] = 1
 128         }
 129         return m
 130 }
 131
 132 // hessrand returns a Hessenberg matrix of order n with random non-zero entries
 133 // in column-major format.
 134 func hessrand(n int, rnd *rand.Rand) []float64 {
 135         h := make([]float64, n*n)
 136         for j := 0; j < n; j++ {
 137                 for i := 0; i <= min(j+1, n-1); i++ {
 138                         h[i+j*n] = rnd.NormFloat64()
 139                 }
 140         }
 141         return h
 142 }
 143
 144 // shiftpairs generates k real and complex conjugate shift pairs. That is, the
 145 // length of sr and si is 2*k.
 146 func shiftpairs(k int, rnd *rand.Rand) (sr, si []float64) {
 147         sr = make([]float64, 2*k)
 148         si = make([]float64, 2*k)
 149         for i := 0; i < len(sr); {
 150                 if rnd.Float64() < 0.5 || i == len(sr)-1 {
 151                         sr[i] = rnd.NormFloat64()
 152                         i++
 153                         continue
 154                 }
 155                 // Generate a complex conjugate pair.
 156                 r := rnd.NormFloat64()
 157                 c := rnd.NormFloat64()
 158                 sr[i] = r
 159                 si[i] = c
 160                 sr[i+1] = r
 161                 si[i+1] = -c
 162                 i += 2
 163         }
 164         return sr, si
 165 }
 166
 167 // rowMajor returns the given r×c column-major matrix a in row-major format.
 168 func rowMajor(r, c int, a []float64) []float64 {
 169         if len(a) != r*c {
 170                 panic("testdata: slice length mismatch")
 171         }
 172         m := make([]float64, len(a))
 173         for i := 0; i < r; i++ {
 174                 for j := 0; j < c; j++ {
 175                         m[i*c+j] = a[i+j*r]
 176                 }
 177         }
 178         return m
 179 }
 180
 181 func min(a, b int) int {
 182         if a < b {
 183                 return a
 184         }
 185         return b
 186 }