1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
// Copyright ©2015 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package testlapack
import (
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
)
type Dgetrier interface {
Dgetrfer
Dgetri(n int, a []float64, lda int, ipiv []int, work []float64, lwork int) bool
}
func DgetriTest(t *testing.T, impl Dgetrier) {
const tol = 1e-13
rnd := rand.New(rand.NewSource(1))
bi := blas64.Implementation()
for _, test := range []struct {
n, lda int
}{
{5, 0},
{5, 8},
{45, 0},
{45, 50},
{63, 70},
{64, 70},
{65, 0},
{65, 70},
{66, 70},
{150, 0},
{150, 250},
} {
n := test.n
lda := test.lda
if lda == 0 {
lda = n
}
// Generate a random well conditioned matrix
perm := rnd.Perm(n)
a := make([]float64, n*lda)
for i := 0; i < n; i++ {
a[i*lda+perm[i]] = 1
}
for i := range a {
a[i] += 0.01 * rnd.Float64()
}
aCopy := make([]float64, len(a))
copy(aCopy, a)
ipiv := make([]int, n)
// Compute LU decomposition.
impl.Dgetrf(n, n, a, lda, ipiv)
// Test with various workspace sizes.
for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
ainv := make([]float64, len(a))
copy(ainv, a)
var lwork int
switch wl {
case minimumWork:
lwork = max(1, n)
case mediumWork:
work := make([]float64, 1)
impl.Dgetri(n, ainv, lda, ipiv, work, -1)
lwork = max(int(work[0])-2*n, n)
case optimumWork:
work := make([]float64, 1)
impl.Dgetri(n, ainv, lda, ipiv, work, -1)
lwork = int(work[0])
}
work := make([]float64, lwork)
// Compute inverse.
ok := impl.Dgetri(n, ainv, lda, ipiv, work, lwork)
if !ok {
t.Errorf("Unexpected singular matrix.")
}
// Check that A(inv) * A = I.
ans := make([]float64, len(ainv))
bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, aCopy, lda, ainv, lda, 0, ans, lda)
// The tolerance is so high because computing matrix inverses is very unstable.
dist := distFromIdentity(n, ans, lda)
if dist > tol {
t.Errorf("|Inv(A) * A - I|_inf = %v is too large. n = %v, lda = %v", dist, n, lda)
}
}
}
}
|
