diff --git a/lapack/testlapack/dgetri.go b/lapack/testlapack/dgetri.go index 4ae08436..bcae50cd 100644 --- a/lapack/testlapack/dgetri.go +++ b/lapack/testlapack/dgetri.go @@ -19,6 +19,7 @@ type Dgetrier interface { } func DgetriTest(t *testing.T, impl Dgetrier) { + const tol = 1e-13 rnd := rand.New(rand.NewSource(1)) bi := blas64.Implementation() for _, test := range []struct { @@ -28,8 +29,11 @@ func DgetriTest(t *testing.T, impl Dgetrier) { {5, 8}, {45, 0}, {45, 50}, + {63, 70}, + {64, 70}, {65, 0}, {65, 70}, + {66, 70}, {150, 0}, {150, 250}, } { @@ -67,8 +71,9 @@ func DgetriTest(t *testing.T, impl Dgetrier) { ans := make([]float64, len(a)) bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, aCopy, lda, a, lda, 0, ans, lda) // The tolerance is so high because computing matrix inverses is very unstable. - if !isIdentity(n, ans, lda, 5e-2) { - t.Errorf("Inv(A) * A != I. n = %v, lda = %v", n, lda) + 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) } } } diff --git a/lapack/testlapack/dlarfg.go b/lapack/testlapack/dlarfg.go index dbed5d4d..1eef8c4a 100644 --- a/lapack/testlapack/dlarfg.go +++ b/lapack/testlapack/dlarfg.go @@ -104,8 +104,9 @@ func DlarfgTest(t *testing.T, impl Dlarfger) { Data: make([]float64, n*n), } blas64.Gemm(blas.Trans, blas.NoTrans, 1, hmat, hmat, 0, eye) - if !isIdentity(n, eye.Data, n, 1e-14) { - t.Errorf("H^T * H is not I %v", eye) + dist := distFromIdentity(n, eye.Data, n) + if dist > 1e-14 { + t.Errorf("H^T * H is not close to I, dist=%v", dist) } xVec := blas64.Vector{ diff --git a/lapack/testlapack/dpotri.go b/lapack/testlapack/dpotri.go index f25d3f4d..7f218037 100644 --- a/lapack/testlapack/dpotri.go +++ b/lapack/testlapack/dpotri.go @@ -97,9 +97,10 @@ func DpotriTest(t *testing.T, impl Dpotrier) { want := make([]float64, n*ldwant) bi.Dsymm(blas.Left, uplo, n, n, 1, aCopy, lda, ainv, ldainv, 0, want, ldwant) - // Check that want is the identity matrix. - if !isIdentity(n, want, ldwant, tol) { - t.Errorf("%v: A * inv(A) != I", prefix) + // Check that want is close to the identity matrix. + dist := distFromIdentity(n, want, ldwant) + if dist > tol { + t.Errorf("%v: |A * inv(A) - I| = %v is too large", prefix, dist) } } } diff --git a/lapack/testlapack/dtrti2.go b/lapack/testlapack/dtrti2.go index 40a4b874..bfaf8a97 100644 --- a/lapack/testlapack/dtrti2.go +++ b/lapack/testlapack/dtrti2.go @@ -147,9 +147,10 @@ func Dtrti2Test(t *testing.T, impl Dtrti2er) { // Compute A^{-1} * A and store the result in ans. ans := make([]float64, len(a)) bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, a, lda, aCopy, lda, 0, ans, lda) - // Check that ans is the identity matrix. - if !isIdentity(n, ans, lda, tol) { - t.Errorf("inv(A) * A != I. Upper = %v, unit = %v, ans = %v", uplo == blas.Upper, diag == blas.Unit, ans) + // Check that ans is close to the identity matrix. + dist := distFromIdentity(n, ans, lda) + if dist > tol { + t.Errorf("|inv(A) * A - I| = %v. Upper = %v, unit = %v, ans = %v", dist, uplo == blas.Upper, diag == blas.Unit, ans) } } } diff --git a/lapack/testlapack/dtrtri.go b/lapack/testlapack/dtrtri.go index 4b0345de..eb8ae648 100644 --- a/lapack/testlapack/dtrtri.go +++ b/lapack/testlapack/dtrtri.go @@ -79,9 +79,10 @@ func DtrtriTest(t *testing.T, impl Dtrtrier) { ans := make([]float64, len(a)) bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, a, lda, aCopy, lda, 0, ans, lda) // Check that ans is the identity matrix. - if !isIdentity(n, ans, lda, tol) { - t.Errorf("inv(A) * A != I. Upper = %v, unit = %v, n = %v, lda = %v", - uplo == blas.Upper, diag == blas.Unit, n, lda) + dist := distFromIdentity(n, ans, lda) + if dist > tol { + t.Errorf("|inv(A) * A - I| = %v is too large. Upper = %v, unit = %v, n = %v, lda = %v", + dist, uplo == blas.Upper, diag == blas.Unit, n, lda) } } } diff --git a/lapack/testlapack/general.go b/lapack/testlapack/general.go index 17ce8ecc..c2373239 100644 --- a/lapack/testlapack/general.go +++ b/lapack/testlapack/general.go @@ -1465,29 +1465,24 @@ func constructGSVPresults(n, p, m, k, l int, a, b blas64.General) (zeroA, zeroB return zeroA, zeroB } -// isIdentity returns whether an n×n matrix A is approximately equal to the -// identity matrix. -func isIdentity(n int, a []float64, lda int, tol float64) bool { +// distFromIdentity returns the L-infinity distance of an n×n matrix A from the +// identity. If A contains NaN elements, distFromIdentity will return +inf. +func distFromIdentity(n int, a []float64, lda int) float64 { + var dist float64 for i := 0; i < n; i++ { for j := 0; j < n; j++ { aij := a[i*lda+j] if math.IsNaN(aij) { - return false + return math.Inf(1) } if i == j { - if math.Abs(aij-1) > tol { - fmt.Println(i, j, aij) - return false - } + dist = math.Max(dist, math.Abs(aij-1)) } else { - if math.Abs(aij) > tol { - fmt.Println(i, j, aij) - return false - } + dist = math.Max(dist, math.Abs(aij)) } } } - return true + return dist } func sameFloat64(a, b float64) bool { diff --git a/lapack/testlapack/matgen_test.go b/lapack/testlapack/matgen_test.go index 0e77ca02..1257487a 100644 --- a/lapack/testlapack/matgen_test.go +++ b/lapack/testlapack/matgen_test.go @@ -32,8 +32,9 @@ func TestDlagsy(t *testing.T) { Dlagsy(n, 0, d, a, lda, rnd, work) // A should be the identity matrix because // A = U * D * U^T = U * I * U^T = U * U^T = I. - if !isIdentity(n, a, lda, tol) { - t.Errorf("Case n=%v,lda=%v: unexpected result", n, lda) + dist := distFromIdentity(n, a, lda) + if dist > tol { + t.Errorf("Case n=%v,lda=%v: |A-I|=%v is too large", n, lda, dist) } } }