lapack/gonum: fix workspace usage in Dgetri
... and extend its test which uncovered these bugs.
MODIFY
lapack/gonum/dgetri.go
lapack/gonum/dgetri.go
@@ -59,14 +59,14 @@ func (impl Implementation) Dgetri(n int, a []float64, lda int, ipiv []int, work
}
nbmin := 2
- ldwork := nb
if 1 < nb && nb < n {
- iws := max(n*ldwork, 1)
+ iws := max(n*nb, 1)
if lwork < iws {
nb = lwork / n
nbmin = max(2, impl.Ilaenv(2, "DGETRI", " ", n, -1, -1, -1))
}
}
+ ldwork := nb
bi := blas64.Implementation()
// Solve the equation inv(A)*L = inv(U) for inv(A).
@@ -76,12 +76,12 @@ func (impl Implementation) Dgetri(n int, a []float64, lda int, ipiv []int, work
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]
+ work[i] = a[i*lda+j]
a[i*lda+j] = 0
}
// 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)
+ bi.Dgemv(blas.NoTrans, n, n-j-1, -1, a[(j+1):], lda, work[(j+1):], 1, 1, a[j:], lda)
}
}
} else {

MODIFY
lapack/testlapack/dgetri.go
lapack/testlapack/dgetri.go
@@ -56,20 +56,38 @@ func DgetriTest(t *testing.T, impl Dgetrier) {
ipiv := make([]int, n)
// Compute LU decomposition.
impl.Dgetrf(n, n, a, lda, ipiv)
- // Compute inverse.
+ // Test with various workspace sizes.
+ for size := range []int{0, 1, 2} {
+ ainv := make([]float64, len(a))
+ copy(ainv, a)
+
+ var lwork int
+ switch size {
+ case 0:
+ // Minimum workspace size.
+ lwork = max(1, n)
+ case 1:
+ // Medium workspace size.
+ work := make([]float64, 1)
+ impl.Dgetri(n, ainv, lda, ipiv, work, -1)
+ lwork = max(int(work[0])-2*n, n)
+ case 2:
+ // Optimum workspace size.
work := make([]float64, 1)
- impl.Dgetri(n, a, lda, ipiv, work, -1)
- work = make([]float64, int(work[0]))
- lwork := len(work)
+ impl.Dgetri(n, ainv, lda, ipiv, work, -1)
+ lwork = int(work[0])
+ }
+ work := make([]float64, lwork)
- ok := impl.Dgetri(n, a, lda, ipiv, work, 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(a))
- bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, aCopy, lda, a, lda, 0, ans, lda)
+ 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 {
@@ -77,3 +95,4 @@ func DgetriTest(t *testing.T, impl Dgetrier) {
}
}
}
+}

lapack/gonum: fix workspace usage in Dgetri

... and extend its test which uncovered these bugs.