lapack/testlapack: change isIdentity to distFromIdentity

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)
 		}
 	}
 }

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{

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)
 					}
 				}
 			}

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)
 				}
 			}
 		}

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)
 				}
 			}
 		}

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 {

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)
 			}
 		}
 	}