vendor/gonum.org/v1/gonum/lapack/testlapack/dpotrf.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         "testing"
   9
  10         "golang.org/x/exp/rand"
  11
  12         "gonum.org/v1/gonum/blas"
  13         "gonum.org/v1/gonum/blas/blas64"
  14         "gonum.org/v1/gonum/floats"
  15 )
  16
  17 type Dpotrfer interface {
  18         Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok bool)
  19 }
  20
  21 func DpotrfTest(t *testing.T, impl Dpotrfer) {
  22         const tol = 1e-13
  23         rnd := rand.New(rand.NewSource(1))
  24         bi := blas64.Implementation()
  25         for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
  26                 for tc, test := range []struct {
  27                         n   int
  28                         lda int
  29                 }{
  30                         {1, 0},
  31                         {2, 0},
  32                         {3, 0},
  33                         {10, 0},
  34                         {30, 0},
  35                         {63, 0},
  36                         {65, 0},
  37                         {127, 0},
  38                         {129, 0},
  39                         {500, 0},
  40                         {1, 10},
  41                         {2, 10},
  42                         {3, 10},
  43                         {10, 20},
  44                         {30, 50},
  45                         {63, 100},
  46                         {65, 100},
  47                         {127, 200},
  48                         {129, 200},
  49                         {500, 600},
  50                 } {
  51                         n := test.n
  52
  53                         // Random diagonal matrix D with positive entries.
  54                         d := make([]float64, n)
  55                         Dlatm1(d, 4, 10000, false, 1, rnd)
  56
  57                         // Construct a positive definite matrix A as
  58                         //  A = U * D * U^T
  59                         // where U is a random orthogonal matrix.
  60                         lda := test.lda
  61                         if lda == 0 {
  62                                 lda = n
  63                         }
  64                         a := make([]float64, n*lda)
  65                         Dlagsy(n, 0, d, a, lda, rnd, make([]float64, 2*n))
  66
  67                         aCopy := make([]float64, len(a))
  68                         copy(aCopy, a)
  69
  70                         ok := impl.Dpotrf(uplo, n, a, lda)
  71                         if !ok {
  72                                 t.Errorf("Case %v: unexpected failure for positive definite matrix", tc)
  73                                 continue
  74                         }
  75
  76                         switch uplo {
  77                         case blas.Upper:
  78                                 for i := 0; i < n; i++ {
  79                                         for j := 0; j < i; j++ {
  80                                                 a[i*lda+j] = 0
  81                                         }
  82                                 }
  83                         case blas.Lower:
  84                                 for i := 0; i < n; i++ {
  85                                         for j := i + 1; j < n; j++ {
  86                                                 a[i*lda+j] = 0
  87                                         }
  88                                 }
  89                         default:
  90                                 panic("bad uplo")
  91                         }
  92
  93                         ans := make([]float64, len(a))
  94                         switch uplo {
  95                         case blas.Upper:
  96                                 // Multiply U^T * U.
  97                                 bi.Dsyrk(uplo, blas.Trans, n, n, 1, a, lda, 0, ans, lda)
  98                         case blas.Lower:
  99                                 // Multiply L * L^T.
 100                                 bi.Dsyrk(uplo, blas.NoTrans, n, n, 1, a, lda, 0, ans, lda)
 101                         }
 102
 103                         match := true
 104                         switch uplo {
 105                         case blas.Upper:
 106                                 for i := 0; i < n; i++ {
 107                                         for j := i; j < n; j++ {
 108                                                 if !floats.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
 109                                                         match = false
 110                                                 }
 111                                         }
 112                                 }
 113                         case blas.Lower:
 114                                 for i := 0; i < n; i++ {
 115                                         for j := 0; j <= i; j++ {
 116                                                 if !floats.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
 117                                                         match = false
 118                                                 }
 119                                         }
 120                                 }
 121                         }
 122                         if !match {
 123                                 t.Errorf("Case %v (uplo=%v,n=%v,lda=%v): unexpected result", tc, uplo, n, lda)
 124                         }
 125
 126                         // Make one element of D negative so that A is not
 127                         // positive definite, and check that Dpotrf fails.
 128                         d[0] *= -1
 129                         Dlagsy(n, 0, d, a, lda, rnd, make([]float64, 2*n))
 130                         ok = impl.Dpotrf(uplo, n, a, lda)
 131                         if ok {
 132                                 t.Errorf("Case %v: unexpected success for not positive definite matrix", tc)
 133                         }
 134                 }
 135         }
 136 }