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.
11 "golang.org/x/exp/rand"
13 "gonum.org/v1/gonum/blas"
14 "gonum.org/v1/gonum/blas/blas64"
17 type Dgetrier interface {
19 Dgetri(n int, a []float64, lda int, ipiv []int, work []float64, lwork int) bool
22 func DgetriTest(t *testing.T, impl Dgetrier) {
23 rnd := rand.New(rand.NewSource(1))
24 bi := blas64.Implementation()
25 for _, test := range []struct {
42 // Generate a random well conditioned matrix
44 a := make([]float64, n*lda)
45 for i := 0; i < n; i++ {
49 a[i] += 0.01 * rnd.Float64()
51 aCopy := make([]float64, len(a))
53 ipiv := make([]int, n)
54 // Compute LU decomposition.
55 impl.Dgetrf(n, n, a, lda, ipiv)
57 work := make([]float64, 1)
58 impl.Dgetri(n, a, lda, ipiv, work, -1)
59 work = make([]float64, int(work[0]))
62 ok := impl.Dgetri(n, a, lda, ipiv, work, lwork)
64 t.Errorf("Unexpected singular matrix.")
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)
71 for i := 0; i < n; i++ {
72 for j := 0; j < n; j++ {
74 // This tolerance is so high because computing matrix inverses
76 if math.Abs(ans[i*lda+j]-1) > 5e-2 {
80 if math.Abs(ans[i*lda+j]) > 5e-2 {
87 t.Errorf("Inv(A) * A != I. n = %v, lda = %v", n, lda)