vendor/gonum.org/v1/gonum/lapack/testlapack/dgetri.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 )
  16
  17 type Dgetrier interface {
  18         Dgetrfer
  19         Dgetri(n int, a []float64, lda int, ipiv []int, work []float64, lwork int) bool
  20 }
  21
  22 func DgetriTest(t *testing.T, impl Dgetrier) {
  23         rnd := rand.New(rand.NewSource(1))
  24         bi := blas64.Implementation()
  25         for _, test := range []struct {
  26                 n, lda int
  27         }{
  28                 {5, 0},
  29                 {5, 8},
  30                 {45, 0},
  31                 {45, 50},
  32                 {65, 0},
  33                 {65, 70},
  34                 {150, 0},
  35                 {150, 250},
  36         } {
  37                 n := test.n
  38                 lda := test.lda
  39                 if lda == 0 {
  40                         lda = n
  41                 }
  42                 // Generate a random well conditioned matrix
  43                 perm := rnd.Perm(n)
  44                 a := make([]float64, n*lda)
  45                 for i := 0; i < n; i++ {
  46                         a[i*lda+perm[i]] = 1
  47                 }
  48                 for i := range a {
  49                         a[i] += 0.01 * rnd.Float64()
  50                 }
  51                 aCopy := make([]float64, len(a))
  52                 copy(aCopy, a)
  53                 ipiv := make([]int, n)
  54                 // Compute LU decomposition.
  55                 impl.Dgetrf(n, n, a, lda, ipiv)
  56                 // Compute inverse.
  57                 work := make([]float64, 1)
  58                 impl.Dgetri(n, a, lda, ipiv, work, -1)
  59                 work = make([]float64, int(work[0]))
  60                 lwork := len(work)
  61
  62                 ok := impl.Dgetri(n, a, lda, ipiv, work, lwork)
  63                 if !ok {
  64                         t.Errorf("Unexpected singular matrix.")
  65                 }
  66
  67                 // Check that A(inv) * A = I.
  68                 ans := make([]float64, len(a))
  69                 bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, aCopy, lda, a, lda, 0, ans, lda)
  70                 isEye := true
  71                 for i := 0; i < n; i++ {
  72                         for j := 0; j < n; j++ {
  73                                 if i == j {
  74                                         // This tolerance is so high because computing matrix inverses
  75                                         // is very unstable.
  76                                         if math.Abs(ans[i*lda+j]-1) > 5e-2 {
  77                                                 isEye = false
  78                                         }
  79                                 } else {
  80                                         if math.Abs(ans[i*lda+j]) > 5e-2 {
  81                                                 isEye = false
  82                                         }
  83                                 }
  84                         }
  85                 }
  86                 if !isEye {
  87                         t.Errorf("Inv(A) * A != I. n = %v, lda = %v", n, lda)
  88                 }
  89         }
  90 }