vendor/gonum.org/v1/gonum/lapack/testlapack/dtrcon.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/floats"
  15         "gonum.org/v1/gonum/lapack"
  16 )
  17
  18 type Dtrconer interface {
  19         Dgeconer
  20         Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64
  21 }
  22
  23 func DtrconTest(t *testing.T, impl Dtrconer) {
  24         rnd := rand.New(rand.NewSource(1))
  25         // Hand crafted tests.
  26         for _, test := range []struct {
  27                 a       []float64
  28                 n       int
  29                 uplo    blas.Uplo
  30                 diag    blas.Diag
  31                 condOne float64
  32                 condInf float64
  33         }{
  34                 {
  35                         a: []float64{
  36                                 8, 5, 6,
  37                                 0, 7, 8,
  38                                 0, 0, 6,
  39                         },
  40                         n:       3,
  41                         uplo:    blas.Upper,
  42                         diag:    blas.Unit,
  43                         condOne: 1.0 / 645,
  44                         condInf: 1.0 / 480,
  45                 },
  46                 {
  47                         a: []float64{
  48                                 8, 5, 6,
  49                                 0, 7, 8,
  50                                 0, 0, 6,
  51                         },
  52                         n:       3,
  53                         uplo:    blas.Upper,
  54                         diag:    blas.NonUnit,
  55                         condOne: 0.137704918032787,
  56                         condInf: 0.157894736842105,
  57                 },
  58                 {
  59                         a: []float64{
  60                                 8, 0, 0,
  61                                 5, 7, 0,
  62                                 6, 8, 6,
  63                         },
  64                         n:       3,
  65                         uplo:    blas.Lower,
  66                         diag:    blas.Unit,
  67                         condOne: 1.0 / 480,
  68                         condInf: 1.0 / 645,
  69                 },
  70                 {
  71                         a: []float64{
  72                                 8, 0, 0,
  73                                 5, 7, 0,
  74                                 6, 8, 6,
  75                         },
  76                         n:       3,
  77                         uplo:    blas.Lower,
  78                         diag:    blas.NonUnit,
  79                         condOne: 0.157894736842105,
  80                         condInf: 0.137704918032787,
  81                 },
  82         } {
  83                 lda := test.n
  84                 work := make([]float64, 3*test.n)
  85                 for i := range work {
  86                         work[i] = rnd.Float64()
  87                 }
  88                 iwork := make([]int, test.n)
  89                 for i := range iwork {
  90                         iwork[i] = int(rnd.Int31())
  91                 }
  92                 aCopy := make([]float64, len(test.a))
  93                 copy(aCopy, test.a)
  94                 condOne := impl.Dtrcon(lapack.MaxColumnSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork)
  95                 if math.Abs(condOne-test.condOne) > 1e-14 {
  96                         t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne)
  97                 }
  98                 if !floats.Equal(aCopy, test.a) {
  99                         t.Errorf("a modified during call")
 100                 }
 101                 condInf := impl.Dtrcon(lapack.MaxRowSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork)
 102                 if math.Abs(condInf-test.condInf) > 1e-14 {
 103                         t.Errorf("Inf norm mismatch. Want %v, got %v.", test.condInf, condInf)
 104                 }
 105                 if !floats.Equal(aCopy, test.a) {
 106                         t.Errorf("a modified during call")
 107                 }
 108         }
 109
 110         // Dtrcon does not match the Dgecon output in many cases. See
 111         // https://github.com/xianyi/OpenBLAS/issues/636
 112         // TODO(btracey): Uncomment this when the mismatch between Dgecon and Dtrcon
 113         // is understood.
 114         /*
 115                 // Randomized tests against Dgecon.
 116                 for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
 117                         for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
 118                                 for _, test := range []struct {
 119                                         n, lda int
 120                                 }{
 121                                         {3, 0},
 122                                         {4, 9},
 123                                 } {
 124                                         for trial := 0; trial < 1; trial++ {
 125                                                 n := test.n
 126                                                 lda := test.lda
 127                                                 if lda == 0 {
 128                                                         lda = n
 129                                                 }
 130                                                 a := make([]float64, n*lda)
 131                                                 if trial == 0 {
 132                                                         for i := range a {
 133                                                                 a[i] = float64(i + 2)
 134                                                         }
 135                                                 } else {
 136                                                         for i := range a {
 137                                                                 a[i] = rnd.NormFloat64()
 138                                                         }
 139                                                 }
 140
 141                                                 aDense := make([]float64, len(a))
 142                                                 if uplo == blas.Upper {
 143                                                         for i := 0; i < n; i++ {
 144                                                                 for j := i; j < n; j++ {
 145                                                                         aDense[i*lda+j] = a[i*lda+j]
 146                                                                 }
 147                                                         }
 148                                                 } else {
 149                                                         for i := 0; i < n; i++ {
 150                                                                 for j := 0; j <= i; j++ {
 151                                                                         aDense[i*lda+j] = a[i*lda+j]
 152                                                                 }
 153                                                         }
 154                                                 }
 155                                                 if diag == blas.Unit {
 156                                                         for i := 0; i < n; i++ {
 157                                                                 aDense[i*lda+i] = 1
 158                                                         }
 159                                                 }
 160
 161                                                 ipiv := make([]int, n)
 162                                                 work := make([]float64, 4*n)
 163                                                 denseOne := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work)
 164                                                 denseInf := impl.Dlange(lapack.MaxRowSum, n, n, aDense, lda, work)
 165
 166                                                 aDenseLU := make([]float64, len(aDense))
 167                                                 copy(aDenseLU, aDense)
 168                                                 impl.Dgetrf(n, n, aDenseLU, lda, ipiv)
 169                                                 iwork := make([]int, n)
 170                                                 want := impl.Dgecon(lapack.MaxColumnSum, n, aDenseLU, lda, denseOne, work, iwork)
 171                                                 got := impl.Dtrcon(lapack.MaxColumnSum, uplo, diag, n, a, lda, work, iwork)
 172                                                 if math.Abs(want-got) > 1e-14 {
 173                                                         t.Errorf("One norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got)
 174                                                 }
 175                                                 want = impl.Dgecon(lapack.MaxRowSum, n, aDenseLU, lda, denseInf, work, iwork)
 176                                                 got = impl.Dtrcon(lapack.MaxRowSum, uplo, diag, n, a, lda, work, iwork)
 177                                                 if math.Abs(want-got) > 1e-14 {
 178                                                         t.Errorf("Inf norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got)
 179                                                 }
 180                                         }
 181                                 }
 182                         }
 183                 }
 184         */
 185 }