vendor/gonum.org/v1/gonum/lapack/testlapack/dtrti2.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/blas/blas64"
  15         "gonum.org/v1/gonum/floats"
  16 )
  17
  18 type Dtrti2er interface {
  19         Dtrti2(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int)
  20 }
  21
  22 func Dtrti2Test(t *testing.T, impl Dtrti2er) {
  23         const tol = 1e-14
  24         for _, test := range []struct {
  25                 a    []float64
  26                 n    int
  27                 uplo blas.Uplo
  28                 diag blas.Diag
  29                 ans  []float64
  30         }{
  31                 {
  32                         a: []float64{
  33                                 2, 3, 4,
  34                                 0, 5, 6,
  35                                 8, 0, 8},
  36                         n:    3,
  37                         uplo: blas.Upper,
  38                         diag: blas.NonUnit,
  39                         ans: []float64{
  40                                 0.5, -0.3, -0.025,
  41                                 0, 0.2, -0.15,
  42                                 8, 0, 0.125,
  43                         },
  44                 },
  45                 {
  46                         a: []float64{
  47                                 5, 3, 4,
  48                                 0, 7, 6,
  49                                 10, 0, 8},
  50                         n:    3,
  51                         uplo: blas.Upper,
  52                         diag: blas.Unit,
  53                         ans: []float64{
  54                                 5, -3, 14,
  55                                 0, 7, -6,
  56                                 10, 0, 8,
  57                         },
  58                 },
  59                 {
  60                         a: []float64{
  61                                 2, 0, 0,
  62                                 3, 5, 0,
  63                                 4, 6, 8},
  64                         n:    3,
  65                         uplo: blas.Lower,
  66                         diag: blas.NonUnit,
  67                         ans: []float64{
  68                                 0.5, 0, 0,
  69                                 -0.3, 0.2, 0,
  70                                 -0.025, -0.15, 0.125,
  71                         },
  72                 },
  73                 {
  74                         a: []float64{
  75                                 1, 0, 0,
  76                                 3, 1, 0,
  77                                 4, 6, 1},
  78                         n:    3,
  79                         uplo: blas.Lower,
  80                         diag: blas.Unit,
  81                         ans: []float64{
  82                                 1, 0, 0,
  83                                 -3, 1, 0,
  84                                 14, -6, 1,
  85                         },
  86                 },
  87         } {
  88                 impl.Dtrti2(test.uplo, test.diag, test.n, test.a, test.n)
  89                 if !floats.EqualApprox(test.ans, test.a, tol) {
  90                         t.Errorf("Matrix inverse mismatch. Want %v, got %v.", test.ans, test.a)
  91                 }
  92         }
  93         rnd := rand.New(rand.NewSource(1))
  94         bi := blas64.Implementation()
  95         for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
  96                 for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
  97                         for _, test := range []struct {
  98                                 n, lda int
  99                         }{
 100                                 {1, 0},
 101                                 {2, 0},
 102                                 {3, 0},
 103                                 {1, 5},
 104                                 {2, 5},
 105                                 {3, 5},
 106                         } {
 107                                 n := test.n
 108                                 lda := test.lda
 109                                 if lda == 0 {
 110                                         lda = n
 111                                 }
 112                                 a := make([]float64, n*lda)
 113                                 for i := range a {
 114                                         a[i] = rnd.Float64()
 115                                 }
 116                                 aCopy := make([]float64, len(a))
 117                                 copy(aCopy, a)
 118                                 impl.Dtrti2(uplo, diag, n, a, lda)
 119                                 if uplo == blas.Upper {
 120                                         for i := 1; i < n; i++ {
 121                                                 for j := 0; j < i; j++ {
 122                                                         aCopy[i*lda+j] = 0
 123                                                         a[i*lda+j] = 0
 124                                                 }
 125                                         }
 126                                 } else {
 127                                         for i := 0; i < n; i++ {
 128                                                 for j := i + 1; j < n; j++ {
 129                                                         aCopy[i*lda+j] = 0
 130                                                         a[i*lda+j] = 0
 131                                                 }
 132                                         }
 133                                 }
 134                                 if diag == blas.Unit {
 135                                         for i := 0; i < n; i++ {
 136                                                 a[i*lda+i] = 1
 137                                                 aCopy[i*lda+i] = 1
 138                                         }
 139                                 }
 140                                 ans := make([]float64, len(a))
 141                                 bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, a, lda, aCopy, lda, 0, ans, lda)
 142                                 iseye := true
 143                                 for i := 0; i < n; i++ {
 144                                         for j := 0; j < n; j++ {
 145                                                 if i == j {
 146                                                         if math.Abs(ans[i*lda+i]-1) > tol {
 147                                                                 iseye = false
 148                                                                 break
 149                                                         }
 150                                                 } else {
 151                                                         if math.Abs(ans[i*lda+j]) > tol {
 152                                                                 iseye = false
 153                                                                 break
 154                                                         }
 155                                                 }
 156                                         }
 157                                 }
 158                                 if !iseye {
 159                                         t.Errorf("inv(A) * A != I. Upper = %v, unit = %v, ans = %v", uplo == blas.Upper, diag == blas.Unit, ans)
 160                                 }
 161                         }
 162                 }
 163         }
 164 }