vendor/gonum.org/v1/gonum/lapack/testlapack/dpocon.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         "log"
   9         "testing"
  10
  11         "golang.org/x/exp/rand"
  12
  13         "gonum.org/v1/gonum/blas"
  14         "gonum.org/v1/gonum/blas/blas64"
  15         "gonum.org/v1/gonum/floats"
  16         "gonum.org/v1/gonum/lapack"
  17 )
  18
  19 type Dpoconer interface {
  20         Dpotrfer
  21         Dgeconer
  22         Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64
  23         Dpocon(uplo blas.Uplo, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
  24 }
  25
  26 func DpoconTest(t *testing.T, impl Dpoconer) {
  27         for _, test := range []struct {
  28                 a    []float64
  29                 n    int
  30                 cond float64
  31                 uplo blas.Uplo
  32         }{
  33                 {
  34                         a: []float64{
  35                                 89, 59, 77,
  36                                 0, 107, 59,
  37                                 0, 0, 89,
  38                         },
  39                         uplo: blas.Upper,
  40                         n:    3,
  41                         cond: 0.050052137643379,
  42                 },
  43                 {
  44                         a: []float64{
  45                                 89, 0, 0,
  46                                 59, 107, 0,
  47                                 77, 59, 89,
  48                         },
  49                         uplo: blas.Lower,
  50                         n:    3,
  51                         cond: 0.050052137643379,
  52                 },
  53                 // Dgecon does not match Dpocon for this case. https://github.com/xianyi/OpenBLAS/issues/664.
  54                 {
  55                         a: []float64{
  56                                 2.9995576045549965, -2.0898894566158663, 3.965560740124006,
  57                                 0, 1.9634729526261008, -2.8681002706874104,
  58                                 0, 0, 5.502416670471008,
  59                         },
  60                         uplo: blas.Upper,
  61                         n:    3,
  62                         cond: 0.024054837369015203,
  63                 },
  64         } {
  65                 n := test.n
  66                 a := make([]float64, len(test.a))
  67                 copy(a, test.a)
  68                 lda := n
  69                 uplo := test.uplo
  70                 work := make([]float64, 3*n)
  71                 anorm := impl.Dlansy(lapack.MaxColumnSum, uplo, n, a, lda, work)
  72                 // Compute cholesky decomposition
  73                 ok := impl.Dpotrf(uplo, n, a, lda)
  74                 if !ok {
  75                         t.Errorf("Bad test, matrix not positive definite")
  76                         continue
  77                 }
  78                 iwork := make([]int, n)
  79                 cond := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork)
  80                 // Error if not the same order, otherwise log the difference.
  81                 if !floats.EqualWithinAbsOrRel(cond, test.cond, 1e0, 1e0) {
  82                         t.Errorf("Cond mismatch. Want %v, got %v.", test.cond, cond)
  83                 } else if !floats.EqualWithinAbsOrRel(cond, test.cond, 1e-14, 1e-14) {
  84                         log.Printf("Dpocon cond mismatch. Want %v, got %v.", test.cond, cond)
  85                 }
  86         }
  87         rnd := rand.New(rand.NewSource(1))
  88         bi := blas64.Implementation()
  89         // Randomized tests compared against Dgecon.
  90         for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
  91                 for _, test := range []struct {
  92                         n, lda int
  93                 }{
  94                         {3, 0},
  95                         {3, 5},
  96                 } {
  97                         for trial := 0; trial < 100; trial++ {
  98                                 n := test.n
  99                                 lda := test.lda
 100                                 if lda == 0 {
 101                                         lda = n
 102                                 }
 103                                 a := make([]float64, n*lda)
 104                                 for i := range a {
 105                                         a[i] = rnd.NormFloat64()
 106                                 }
 107
 108                                 // Multiply a by itself to make it symmetric positive definite.
 109                                 aCopy := make([]float64, len(a))
 110                                 copy(aCopy, a)
 111                                 bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, aCopy, lda, aCopy, lda, 0, a, lda)
 112
 113                                 aDat := make([]float64, len(aCopy))
 114                                 copy(aDat, a)
 115
 116                                 aDense := make([]float64, len(a))
 117                                 if uplo == blas.Upper {
 118                                         for i := 0; i < n; i++ {
 119                                                 for j := i; j < n; j++ {
 120                                                         v := a[i*lda+j]
 121                                                         aDense[i*lda+j] = v
 122                                                         aDense[j*lda+i] = v
 123                                                 }
 124                                         }
 125                                 } else {
 126                                         for i := 0; i < n; i++ {
 127                                                 for j := 0; j <= i; j++ {
 128                                                         v := a[i*lda+j]
 129                                                         aDense[i*lda+j] = v
 130                                                         aDense[j*lda+i] = v
 131                                                 }
 132                                         }
 133                                 }
 134                                 work := make([]float64, 4*n)
 135                                 iwork := make([]int, n)
 136
 137                                 anorm := impl.Dlansy(lapack.MaxColumnSum, uplo, n, a, lda, work)
 138                                 ok := impl.Dpotrf(uplo, n, a, lda)
 139                                 if !ok {
 140                                         t.Errorf("Bad test, matrix not positive definite")
 141                                         continue
 142                                 }
 143                                 got := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork)
 144
 145                                 denseNorm := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work)
 146                                 ipiv := make([]int, n)
 147                                 impl.Dgetrf(n, n, aDense, lda, ipiv)
 148                                 want := impl.Dgecon(lapack.MaxColumnSum, n, aDense, lda, denseNorm, work, iwork)
 149                                 // Error if not the same order, otherwise log the difference.
 150                                 if !floats.EqualWithinAbsOrRel(want, got, 1e0, 1e0) {
 151                                         t.Errorf("Dpocon and Dgecon mismatch. Dpocon %v, Dgecon %v.", got, want)
 152                                 } else if !floats.EqualWithinAbsOrRel(want, got, 1e-14, 1e-14) {
 153                                         log.Printf("Dpocon and Dgecon mismatch. Dpocon %v, Dgecon %v.", got, want)
 154                                 }
 155                         }
 156                 }
 157         }
 158 }