mat: fix Dense.Inverse
The condition number of A was computed by Gecon not from the LU
factorization of A but rather from A^{-1} due to the preceding call to
Getri.
Fixed by reordering the various lapack calls and by allocating
an additional integer work slice.
MODIFY
mat/dense_arithmetic.go
mat/dense_arithmetic.go
@@ -237,28 +237,35 @@ func (m *Dense) Inverse(a Matrix) error {
default:
m.Copy(a)
}
+ // Compute the norm of A.
+ work := getFloats(4*r, false) // Length must be at least 4*r for Gecon.
+ norm := lapack64.Lange(CondNorm, m.mat, work)
+ // Compute the LU factorization of A.
ipiv := getInts(r, false)
defer putInts(ipiv)
ok := lapack64.Getrf(m.mat, ipiv)
if !ok {
+ // A is exactly singular.
return Condition(math.Inf(1))
}
- work := getFloats(4*r, false) // must be at least 4*r for cond.
+ // Compute the condition number of A using the LU factorization.
+ iwork := getInts(r, false)
+ defer putInts(iwork)
+ rcond := lapack64.Gecon(CondNorm, m.mat, norm, work, iwork)
+ // Compute A^{-1} from the LU factorization regardless of the value of rcond.
lapack64.Getri(m.mat, ipiv, work, -1)
- if int(work[0]) > 4*r {
+ if int(work[0]) > len(work) {
l := int(work[0])
putFloats(work)
work = getFloats(l, false)
- } else {
- work = work[:4*r]
}
defer putFloats(work)
- lapack64.Getri(m.mat, ipiv, work, len(work))
- norm := lapack64.Lange(CondNorm, m.mat, work)
- rcond := lapack64.Gecon(CondNorm, m.mat, norm, work, ipiv) // reuse ipiv
- if rcond == 0 {
+ ok = lapack64.Getri(m.mat, ipiv, work, len(work))
+ if !ok || rcond == 0 {
+ // A is exactly singular.
return Condition(math.Inf(1))
}
+ // Check whether A is singular for computational purposes.
cond := 1 / rcond
if cond > ConditionTolerance {
return Condition(cond)

MODIFY
mat/dense_test.go
mat/dense_test.go
@@ -2079,6 +2079,60 @@ func TestDenseInverse(t *testing.T) {
fmt.Println(Formatted(&m))
}
}
+
+ // Randomized tests
+ const tol = 1e-16
+ rnd := rand.New(rand.NewSource(1))
+ for _, recvSameAsA := range []bool{false, true} {
+ for _, trans := range []bool{false, true} {
+ if trans && recvSameAsA {
+ // Transposed argument cannot be a receiver.
+ continue
+ }
+ for _, n := range []int{1, 2, 3, 5, 10, 50, 100} {
+ name := fmt.Sprintf("n=%d,recvSameAsA=%v,trans=%v", n, recvSameAsA, trans)
+
+ // Generate the contents of a random dense matrix.
+ data := make([]float64, n*n)
+ for i := range data {
+ data[i] = rnd.NormFloat64()
+ }
+
+ var a Matrix
+ var aInv Dense
+ if recvSameAsA {
+ aInv = *NewDense(n, n, data)
+ a = &aInv
+ } else {
+ if trans {
+ a = NewDense(n, n, data).T()
+ } else {
+ a = asBasicMatrix(NewDense(n, n, data))
+ }
+
+ }
+ var aOrig Dense
+ aOrig.CloneFrom(a)
+
+ err := aInv.Inverse(a)
+ if err != nil {
+ t.Errorf("%v: unexpected failure of Inverse, %v", name, err)
+ continue
+ }
+
+ // Compute the residual |I - A^{-1} * A| / (N * |A| * |A^{-1}).
+ var aaInv Dense
+ aaInv.Mul(&aInv, &aOrig)
+ for i := 0; i < n; i++ {
+ aaInv.Set(i, i, 1-aaInv.At(i, i))
+ }
+ resid := Norm(&aaInv, 1) / (float64(n) * Norm(&aOrig, 1) * Norm(&aInv, 1))
+ if resid > tol {
+ t.Errorf("%v: A*A^{-1} is not identity, resid=%v,want<=%v", name, resid, tol)
+ }
+ }
+ }
+ }
}
var (

mat: fix Dense.Inverse

The condition number of A was computed by Gecon not from the LU
factorization of A but rather from A^{-1} due to the preceding call to
Getri.

Fixed by reordering the various lapack calls and by allocating
an additional integer work slice.