From bffde37f95a4990c808186401eee72cc79b46a70 Mon Sep 17 00:00:00 2001 From: btracey Date: Wed, 9 Sep 2015 23:00:50 -0600 Subject: [PATCH 1/3] Add Dpocon and tests --- cgo/lapack.go | 23 ++++++++ cgo/lapack_test.go | 4 ++ native/dpocon.go | 72 +++++++++++++++++++++++ native/lapack_test.go | 4 ++ testlapack/dpocon.go | 131 ++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 234 insertions(+) create mode 100644 native/dpocon.go create mode 100644 testlapack/dpocon.go diff --git a/cgo/lapack.go b/cgo/lapack.go index 0d91b20..afd01aa 100644 --- a/cgo/lapack.go +++ b/cgo/lapack.go @@ -524,6 +524,29 @@ func (impl Implementation) Dormqr(side blas.Side, trans blas.Transpose, m, n, k clapack.Dormqr(side, trans, m, n, k, a, lda, tau, c, ldc) } +// Dtrcon estimates the reciprocal of the condition number of a positive-definite +// matrix A given the Cholesky decmposition of A. The condition number computed +// is based on the 1-norm and the ∞-norm. +// +// anorm is the 1-norm and the ∞-norm of the original matrix A. +// +// work is a temporary data slice of length at least 3*n and Dpocon will panic otherwise. +// +// iwork is a temporary data slice of length at least n and Dpocon will panic otherwise. +func (impl Implementation) Dpocon(uplo blas.Uplo, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 { + checkMatrix(n, n, a, lda) + if uplo != blas.Upper && uplo != blas.Lower { + panic(badUplo) + } + if len(work) < 3*n { + panic(badWork) + } + if len(iwork) < n { + panic(badWork) + } + clapack.Dpocon(uplo, n, a, lda, anorm, rcond) +} + // Dtrcon estimates the reciprocal of the condition number of a triangular matrix A. // The condition number computed may be based on the 1-norm or the ∞-norm. // diff --git a/cgo/lapack_test.go b/cgo/lapack_test.go index b4e22c3..7167616 100644 --- a/cgo/lapack_test.go +++ b/cgo/lapack_test.go @@ -87,6 +87,10 @@ func TestDormlq(t *testing.T) { } */ +func TestDpocon(t *testing.T) { + testlapack.DpoconTest(t, impl) +} + func TestDtrcon(t *testing.T) { testlapack.DtrconTest(t, impl) } diff --git a/native/dpocon.go b/native/dpocon.go new file mode 100644 index 0000000..d1ebb78 --- /dev/null +++ b/native/dpocon.go @@ -0,0 +1,72 @@ +package native + +import ( + "math" + + "github.com/gonum/blas" + "github.com/gonum/blas/blas64" +) + +// Dtrcon estimates the reciprocal of the condition number of a positive-definite +// matrix A given the Cholesky decmposition of A. The condition number computed +// is based on the 1-norm and the ∞-norm. +// +// anorm is the 1-norm and the ∞-norm of the original matrix A. +// +// work is a temporary data slice of length at least 3*n and Dpocon will panic otherwise. +// +// iwork is a temporary data slice of length at least n and Dpocon will panic otherwise. +func (impl Implementation) Dpocon(uplo blas.Uplo, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 { + checkMatrix(n, n, a, lda) + if uplo != blas.Upper && uplo != blas.Lower { + panic(badUplo) + } + if len(work) < 3*n { + panic(badWork) + } + if len(iwork) < n { + panic(badWork) + } + var rcond float64 + if n == 0 { + return 1 + } + if anorm == 0 { + return rcond + } + + bi := blas64.Implementation() + var ainvnm float64 + smlnum := dlamchS + upper := uplo == blas.Upper + var kase int + var normin bool + isave := new([3]int) + var sl, su float64 + for { + ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave) + if kase == 0 { + if ainvnm != 0 { + rcond = (1 / ainvnm) / anorm + } + return rcond + } + if upper { + sl = impl.Dlatrs(blas.Upper, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[2*n:]) + normin = true + su = impl.Dlatrs(blas.Upper, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[2*n:]) + } else { + sl = impl.Dlatrs(blas.Lower, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[2*n:]) + normin = true + su = impl.Dlatrs(blas.Lower, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[2*n:]) + } + scale := sl * su + if scale != 1 { + ix := bi.Idamax(n, work, 1) + if scale == 0 || scale < math.Abs(work[ix])*smlnum { + return rcond + } + impl.Drscl(n, scale, work, 1) + } + } +} diff --git a/native/lapack_test.go b/native/lapack_test.go index cadc6ef..9c04df5 100644 --- a/native/lapack_test.go +++ b/native/lapack_test.go @@ -92,6 +92,10 @@ func TestDorm2r(t *testing.T) { testlapack.Dorm2rTest(t, impl) } +func TestDpocon(t *testing.T) { + testlapack.DpoconTest(t, impl) +} + func TestDpotf2(t *testing.T) { testlapack.Dpotf2Test(t, impl) } diff --git a/testlapack/dpocon.go b/testlapack/dpocon.go new file mode 100644 index 0000000..2e5a026 --- /dev/null +++ b/testlapack/dpocon.go @@ -0,0 +1,131 @@ +package testlapack + +import ( + "math" + "math/rand" + "testing" + + "github.com/gonum/blas" + "github.com/gonum/blas/blas64" + "github.com/gonum/lapack" +) + +type Dpoconer interface { + Dpotrfer + Dgeconer + Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64 + Dpocon(uplo blas.Uplo, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 +} + +func DpoconTest(t *testing.T, impl Dpoconer) { + for _, test := range []struct { + a []float64 + n int + cond float64 + uplo blas.Uplo + }{ + { + a: []float64{ + 89, 59, 77, + 0, 107, 59, + 0, 0, 89, + }, + uplo: blas.Upper, + n: 3, + cond: 0.050052137643379, + }, + { + a: []float64{ + 89, 0, 0, + 59, 107, 0, + 77, 59, 89, + }, + uplo: blas.Lower, + n: 3, + cond: 0.050052137643379, + }, + } { + n := test.n + a := make([]float64, len(test.a)) + copy(a, test.a) + lda := n + uplo := test.uplo + work := make([]float64, 3*n) + anorm := impl.Dlansy(lapack.MaxColumnSum, uplo, n, a, lda, work) + // Compute cholesky decomposition + ok := impl.Dpotrf(uplo, n, a, lda) + if !ok { + t.Errorf("Bad test, matrix not positive definite") + continue + } + iwork := make([]int, n) + cond := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork) + if math.Abs(cond-test.cond) > 1e-14 { + t.Errorf("Cond mismatch. Want %v, got %v.", test.cond, cond) + } + } + bi := blas64.Implementation() + // Randomized tests compared against Dgecon. + for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} { + for _, test := range []struct { + n, lda int + }{ + {3, 0}, + {3, 5}, + } { + for trial := 0; trial < 100; trial++ { + n := test.n + lda := test.lda + if lda == 0 { + lda = n + } + a := make([]float64, n*lda) + for i := range a { + a[i] = rand.NormFloat64() + } + + // Multiply a by itself to make it symmetric positive definite. + aCopy := make([]float64, len(a)) + copy(aCopy, a) + bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, aCopy, lda, aCopy, lda, 0, a, lda) + + aDense := make([]float64, len(a)) + if uplo == blas.Upper { + for i := 0; i < n; i++ { + for j := i; j < n; j++ { + v := a[i*lda+j] + aDense[i*lda+j] = v + aDense[j*lda+i] = v + } + } + } else { + for i := 0; i < n; i++ { + for j := 0; j <= i; j++ { + v := a[i*lda+j] + aDense[i*lda+j] = v + aDense[j*lda+i] = v + } + } + } + work := make([]float64, 4*n) + iwork := make([]int, n) + + anorm := impl.Dlansy(lapack.MaxColumnSum, uplo, n, a, lda, work) + ok := impl.Dpotrf(uplo, n, a, lda) + if !ok { + t.Errorf("Bad test, matrix not positive definite") + continue + } + got := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork) + + denseNorm := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work) + ipiv := make([]int, n) + impl.Dgetrf(n, n, aDense, lda, ipiv) + want := impl.Dgecon(lapack.MaxColumnSum, n, aDense, lda, denseNorm, work, iwork) + if math.Abs(got-want) > 1e-14 { + t.Errorf("Cond mismatch. Want %v, got %v.", want, got) + } + } + } + } +} From df3af3a123cd3811f07971c4142164e59aeb84be Mon Sep 17 00:00:00 2001 From: btracey Date: Wed, 9 Sep 2015 23:10:10 -0600 Subject: [PATCH 2/3] Fixed missing return --- cgo/lapack.go | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cgo/lapack.go b/cgo/lapack.go index afd01aa..ad34cbd 100644 --- a/cgo/lapack.go +++ b/cgo/lapack.go @@ -544,7 +544,9 @@ func (impl Implementation) Dpocon(uplo blas.Uplo, n int, a []float64, lda int, a if len(iwork) < n { panic(badWork) } + rcond := make([]float64, 1) clapack.Dpocon(uplo, n, a, lda, anorm, rcond) + return rcond[0] } // Dtrcon estimates the reciprocal of the condition number of a triangular matrix A. From 7d242475878e01be83c0309ecbdeb702d0dd5256 Mon Sep 17 00:00:00 2001 From: btracey Date: Wed, 9 Sep 2015 23:58:49 -0600 Subject: [PATCH 3/3] Disable cgo test --- cgo/lapack_test.go | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cgo/lapack_test.go b/cgo/lapack_test.go index 7167616..2bef153 100644 --- a/cgo/lapack_test.go +++ b/cgo/lapack_test.go @@ -87,9 +87,12 @@ func TestDormlq(t *testing.T) { } */ +/* +// Test disabled because of memory crashes. https://github.com/xianyi/OpenBLAS/issues/639. func TestDpocon(t *testing.T) { testlapack.DpoconTest(t, impl) } +*/ func TestDtrcon(t *testing.T) { testlapack.DtrconTest(t, impl)