lapack/gonum: fix various bugs in Dgetri

lapack/gonum/dgetri.go

@@ -22,71 +22,95 @@ import (
 // by the temporary space available. If lwork == -1, instead of performing Dgetri,
 // the optimal work length will be stored into work[0].
 func (impl Implementation) Dgetri(n int, a []float64, lda int, ipiv []int, work []float64, lwork int) (ok bool) {
-	checkMatrix(n, n, a, lda)
-	if len(ipiv) < n {
-		panic(badIpiv)
+	switch {
+	case n < 0:
+		panic("lapack: has negative number of columns")
+	case lda < max(1, n):
+		panic("lapack: stride less than number of columns")
+	case lwork < max(1, n) && lwork != -1:
+		panic(badWork)
+	case len(work) < max(1, lwork):
+		panic(shortWork)
 	}
-	nb := impl.Ilaenv(1, "DGETRI", " ", n, -1, -1, -1)
-	if lwork == -1 {
-		work[0] = float64(n * nb)
+
+	if n == 0 {
+		work[0] = 1
 		return true
 	}
-	if lwork < n {
-		panic(badWork)
-	}
-	if len(work) < lwork {
-		panic(badWork)
+
+	switch {
+	case len(a) < (n-1)*lda+n:
+		panic("lapack: insufficient matrix slice length")
+	case len(ipiv) < n:
+		panic(badIpiv)
 	}
-	if n == 0 {
+
+	nb := impl.Ilaenv(1, "DGETRI", " ", n, -1, -1, -1)
+	lworkopt := float64(n * nb)
+	if lwork == -1 {
+		work[0] = lworkopt
 		return true
 	}
+
+	// Form inv(U).
 	ok = impl.Dtrtri(blas.Upper, blas.NonUnit, n, a, lda)
 	if !ok {
 		return false
 	}
+
 	nbmin := 2
 	ldwork := nb
-	if nb > 1 && nb < n {
-		iws := max(ldwork*n, 1)
+	if 1 < nb && nb < n {
+		iws := max(n*ldwork, 1)
 		if lwork < iws {
-			nb = lwork / ldwork
+			nb = lwork / n
 			nbmin = max(2, impl.Ilaenv(2, "DGETRI", " ", n, -1, -1, -1))
 		}
 	}
+
 	bi := blas64.Implementation()
+	// Solve the equation inv(A)*L = inv(U) for inv(A).
 	// TODO(btracey): Replace this with a more row-major oriented algorithm.
-	if nb < nbmin || nb >= n {
+	if nb < nbmin || n <= nb {
 		// Unblocked code.
 		for j := n - 1; j >= 0; j-- {
 			for i := j + 1; i < n; i++ {
+				// Copy current column of L to work and replace with zeros.
 				work[i*ldwork] = a[i*lda+j]
 				a[i*lda+j] = 0
 			}
-			if j < n {
+			// Compute current column of inv(A).
+			if j < n-1 {
 				bi.Dgemv(blas.NoTrans, n, n-j-1, -1, a[(j+1):], lda, work[(j+1)*ldwork:], ldwork, 1, a[j:], lda)
 			}
 		}
 	} else {
+		// Blocked code.
 		nn := ((n - 1) / nb) * nb
 		for j := nn; j >= 0; j -= nb {
 			jb := min(nb, n-j)
-			for jj := j; jj < j+jb-1; jj++ {
+			// Copy current block column of L to work and replace
+			// with zeros.
+			for jj := j; jj < j+jb; jj++ {
 				for i := jj + 1; i < n; i++ {
 					work[i*ldwork+(jj-j)] = a[i*lda+jj]
 					a[i*lda+jj] = 0
 				}
 			}
+			// Compute current block column of inv(A).
 			if j+jb < n {
 				bi.Dgemm(blas.NoTrans, blas.NoTrans, n, jb, n-j-jb, -1, a[(j+jb):], lda, work[(j+jb)*ldwork:], ldwork, 1, a[j:], lda)
-				bi.Dtrsm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, n, jb, 1, work[j*ldwork:], ldwork, a[j:], lda)
 			}
+			bi.Dtrsm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, n, jb, 1, work[j*ldwork:], ldwork, a[j:], lda)
 		}
 	}
+	// Apply column interchanges.
 	for j := n - 2; j >= 0; j-- {
 		jp := ipiv[j]
 		if jp != j {
 			bi.Dswap(n, a[j:], lda, a[jp:], lda)
 		}
 	}
+	work[0] = lworkopt
 	return true
 }

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

@@ -12,13 +12,15 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/floats"
 )
 
 type Dlarfger interface {
 	Dlarfg(n int, alpha float64, x []float64, incX int) (beta, tau float64)
 }
 
 func DlarfgTest(t *testing.T, impl Dlarfger) {
+	const tol = 1e-14
 	rnd := rand.New(rand.NewSource(1))
 	for i, test := range []struct {
 		alpha float64
@@ -104,8 +106,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 > tol {
+			t.Errorf("H^T * H is not close to I, dist=%v", dist)
 		}
 
 		xVec := blas64.Vector{
@@ -121,14 +124,11 @@ func DlarfgTest(t *testing.T, impl Dlarfger) {
 			Data: ans,
 		}
 		blas64.Gemv(blas.NoTrans, 1, hmat, xVec, 0, ansVec)
-		if math.Abs(ans[0]-beta) > 1e-14 {
+		if math.Abs(ans[0]-beta) > tol {
 			t.Errorf("Case %v, beta mismatch. Want %v, got %v", i, ans[0], beta)
 		}
-		for i := 1; i < n; i++ {
-			if math.Abs(ans[i]) > 1e-14 {
-				t.Errorf("Case %v, nonzero answer %v", i, ans)
-				break
-			}
+		if floats.Norm(ans[1:n], math.Inf(1)) > tol {
+			t.Errorf("Case %v, nonzero answer %v", i, ans[1:n])
 		}
 	}
 }

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