mat: factor condition code into factorising types
MODIFY
mat/lu.go
mat/lu.go
@@ -10,6 +10,7 @@ import (
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
"gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/lapack"
"gonum.org/v1/gonum/lapack/lapack64"
)
@@ -22,15 +23,15 @@ type LU struct {
cond float64
}
-// updateCond updates the stored condition number of the matrix. Norm is the
-// norm of the original matrix. If norm is negative it will be estimated.
-func (lu *LU) updateCond(norm float64) {
+// updateCond updates the stored condition number of the matrix. anorm is the
+// norm of the original matrix. If anorm is negative it will be estimated.
+func (lu *LU) updateCond(anorm float64, norm lapack.MatrixNorm) {
n := lu.lu.mat.Cols
work := getFloats(4*n, false)
defer putFloats(work)
iwork := getInts(n, false)
defer putInts(iwork)
- if norm < 0 {
+ if anorm < 0 {
// This is an approximation. By the definition of a norm, ||AB|| <= ||A|| ||B||.
// The condition number is ||A|| || A^-1||, so this will underestimate
// the condition number somewhat.
@@ -38,11 +39,11 @@ func (lu *LU) updateCond(norm float64) {
// update possibilities, e.g. RankOne.
u := lu.lu.asTriDense(n, blas.NonUnit, blas.Upper)
l := lu.lu.asTriDense(n, blas.Unit, blas.Lower)
- unorm := lapack64.Lantr(CondNorm, u.mat, work)
- lnorm := lapack64.Lantr(CondNorm, l.mat, work)
- norm = unorm * lnorm
+ unorm := lapack64.Lantr(norm, u.mat, work)
+ lnorm := lapack64.Lantr(norm, l.mat, work)
+ anorm = unorm * lnorm
}
- v := lapack64.Gecon(CondNorm, lu.lu.mat, norm, work, iwork)
+ v := lapack64.Gecon(norm, lu.lu.mat, anorm, work, iwork)
lu.cond = 1 / v
}
@@ -54,6 +55,10 @@ func (lu *LU) updateCond(norm float64) {
// factors can be extracted from the factorization using the Permutation method
// on Dense, and the LU LTo and UTo methods.
func (lu *LU) Factorize(a Matrix) {
+ lu.factorize(a, CondNorm)
+}
+
+func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) {
r, c := a.Dims()
if r != c {
panic(ErrSquare)
@@ -70,10 +75,10 @@ func (lu *LU) Factorize(a Matrix) {
}
lu.pivot = lu.pivot[:r]
work := getFloats(r, false)
- anorm := lapack64.Lange(CondNorm, lu.lu.mat, work)
+ anorm := lapack64.Lange(norm, lu.lu.mat, work)
putFloats(work)
lapack64.Getrf(lu.lu.mat, lu.pivot)
- lu.updateCond(anorm)
+ lu.updateCond(anorm, norm)
}
// Cond returns the condition number for the factorized matrix.
@@ -215,7 +220,7 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y *VecDense) {
lum.Data[j*lum.Stride+i] += gamma * tmp
}
}
- lu.updateCond(-1)
+ lu.updateCond(-1, CondNorm)
}
// LTo extracts the lower triangular matrix from an LU factorization.

MODIFY
mat/matrix.go
mat/matrix.go
@@ -302,63 +302,23 @@ func Cond(a Matrix, norm float64) float64 {
case math.Inf(1):
lnorm = lapack.MaxRowSum
}
+
if m == n {
// Use the LU decomposition to compute the condition number.
- tmp := getWorkspace(m, n, false)
- tmp.Copy(a)
- work := getFloats(4*n, false)
- aNorm := lapack64.Lange(lnorm, tmp.mat, work)
- pivot := getInts(m, false)
- lapack64.Getrf(tmp.mat, pivot)
- iwork := make([]int, n)
- v := lapack64.Gecon(lnorm, tmp.mat, aNorm, work, iwork)
- putWorkspace(tmp)
- putFloats(work)
- putInts(pivot)
- return 1 / v
+ var lu LU
+ lu.factorize(a, lnorm)
+ return lu.Cond()
}
if m > n {
// Use the QR factorization to compute the condition number.
- tmp := getWorkspace(m, n, false)
- tmp.Copy(a)
- work := getFloats(3*n, false)
- tau := getFloats(min(m, n), false)
- lapack64.Geqrf(tmp.mat, tau, work, -1)
- if l := int(work[0]); l > len(work) {
- putFloats(work)
- work = getFloats(l, false)
- }
- lapack64.Geqrf(tmp.mat, tau, work, len(work))
-
- iwork := getInts(n, false)
- r := tmp.asTriDense(n, blas.NonUnit, blas.Upper)
- v := lapack64.Trcon(lnorm, r.mat, work, iwork)
- putWorkspace(tmp)
- putFloats(work)
- putFloats(tau)
- putInts(iwork)
- return 1 / v
+ var qr QR
+ qr.factorize(a, lnorm)
+ return qr.Cond()
}
// Use the LQ factorization to compute the condition number.
- tmp := getWorkspace(m, n, false)
- tmp.Copy(a)
- work := getFloats(3*m, false)
- tau := getFloats(min(m, n), false)
- lapack64.Gelqf(tmp.mat, tau, work, -1)
- if l := int(work[0]); l > len(work) {
- putFloats(work)
- work = getFloats(l, false)
- }
- lapack64.Gelqf(tmp.mat, tau, work, len(work))
-
- iwork := getInts(m, false)
- l := tmp.asTriDense(m, blas.NonUnit, blas.Lower)
- v := lapack64.Trcon(lnorm, l.mat, work, iwork)
- putWorkspace(tmp)
- putFloats(work)
- putFloats(tau)
- putInts(iwork)
- return 1 / v
+ var lq LQ
+ lq.factorize(a, lnorm)
+ return lq.Cond()
}
// Det returns the determinant of the matrix a. In many expressions using LogDet

MODIFY
mat/matrix_test.go
mat/matrix_test.go
@@ -237,18 +237,28 @@ func TestCond(t *testing.T) {
condInf: 31.153846153846136,
},
} {
+ orig := DenseCopyOf(test.a)
condOne := Cond(test.a, 1)
if !floats.EqualWithinAbsOrRel(test.condOne, condOne, 1e-13, 1e-13) {
- t.Errorf("Case %d: inf norm mismatch. Want %v, got %v", i, test.condOne, condOne)
+ t.Errorf("Case %d: one norm mismatch. Want %v, got %v", i, test.condOne, condOne)
+ }
+ if !Equal(test.a, orig) {
+ t.Errorf("Case %d: unexpected mutation of input matrix for one norm. Want %v, got %v", i, orig, test.a)
}
condTwo := Cond(test.a, 2)
if !floats.EqualWithinAbsOrRel(test.condTwo, condTwo, 1e-13, 1e-13) {
- t.Errorf("Case %d: inf norm mismatch. Want %v, got %v", i, test.condTwo, condTwo)
+ t.Errorf("Case %d: two norm mismatch. Want %v, got %v", i, test.condTwo, condTwo)
+ }
+ if !Equal(test.a, orig) {
+ t.Errorf("Case %d: unexpected mutation of input matrix for two norm. Want %v, got %v", i, orig, test.a)
}
condInf := Cond(test.a, math.Inf(1))
if !floats.EqualWithinAbsOrRel(test.condInf, condInf, 1e-13, 1e-13) {
t.Errorf("Case %d: inf norm mismatch. Want %v, got %v", i, test.condInf, condInf)
}
+ if !Equal(test.a, orig) {
+ t.Errorf("Case %d: unexpected mutation of input matrix for inf norm. Want %v, got %v", i, orig, test.a)
+ }
}
for _, test := range []struct {

mat: factor condition code into factorising types