From 7476e69309f1beef9f2a649a2c1974229150000a Mon Sep 17 00:00:00 2001 From: btracey Date: Wed, 9 Sep 2015 21:21:48 -0600 Subject: [PATCH] Add Dtrcon and test --- cgo/lapack.go | 27 +++++++ cgo/lapack_test.go | 4 + native/dtrcon.go | 78 ++++++++++++++++++ native/lapack_test.go | 4 + testlapack/dtrcon.go | 179 ++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 292 insertions(+) create mode 100644 native/dtrcon.go create mode 100644 testlapack/dtrcon.go diff --git a/cgo/lapack.go b/cgo/lapack.go index 63f17101..0d91b208 100644 --- a/cgo/lapack.go +++ b/cgo/lapack.go @@ -524,6 +524,33 @@ 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 triangular matrix A. +// The condition number computed may be based on the 1-norm or the ∞-norm. +// +// work is a temporary data slice of length at least 3*n and Dtrcon will panic otherwise. +// +// iwork is a temporary data slice of length at least n and Dtrcon will panic otherwise. +func (impl Implementation) Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64 { + if norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum { + panic(badNorm) + } + if uplo != blas.Upper && uplo != blas.Lower { + panic(badUplo) + } + if diag != blas.NonUnit && diag != blas.Unit { + panic(badDiag) + } + if len(work) < 3*n { + panic(badWork) + } + if len(iwork) < n { + panic(badWork) + } + rcond := []float64{0} + clapack.Dtrcon(byte(norm), uplo, diag, n, a, lda, rcond) + return rcond[0] +} + // Dtrtrs solves a triangular system of the form A * X = B or A^T * X = B. Dtrtrs // returns whether the solve completed successfully. If A is singular, no solve is performed. func (impl Implementation) Dtrtrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, a []float64, lda int, b []float64, ldb int) (ok bool) { diff --git a/cgo/lapack_test.go b/cgo/lapack_test.go index 14ec3beb..bb667d45 100644 --- a/cgo/lapack_test.go +++ b/cgo/lapack_test.go @@ -89,3 +89,7 @@ func TestDormlq(t *testing.T) { testlapack.Dorml2Test(t, blockedTranslate{impl}) } */ + +func TestDtrcon(t *testing.T) { + testlapack.DtrconTest(t, impl) +} diff --git a/native/dtrcon.go b/native/dtrcon.go new file mode 100644 index 00000000..640fe550 --- /dev/null +++ b/native/dtrcon.go @@ -0,0 +1,78 @@ +package native + +import ( + "math" + + "github.com/gonum/blas" + "github.com/gonum/blas/blas64" + "github.com/gonum/lapack" +) + +// 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. +// +// work is a temporary data slice of length at least 3*n and Dtrcon will panic otherwise. +// +// iwork is a temporary data slice of length at least n and Dtrcon will panic otherwise. +func (impl Implementation) Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64 { + if norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum { + panic(badNorm) + } + if uplo != blas.Upper && uplo != blas.Lower { + panic(badUplo) + } + if diag != blas.NonUnit && diag != blas.Unit { + panic(badDiag) + } + if len(work) < 3*n { + panic(badWork) + } + if len(iwork) < n { + panic(badWork) + } + if n == 0 { + return 1 + } + bi := blas64.Implementation() + + var rcond float64 + smlnum := dlamchS * float64(n) + + anorm := impl.Dlantr(norm, uplo, diag, n, n, a, lda, work) + + if anorm <= 0 { + return rcond + } + var ainvnm float64 + var normin bool + kase1 := 2 + if norm == lapack.MaxColumnSum { + kase1 = 1 + } + var kase int + isave := new([3]int) + var scale float64 + for { + ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave) + if kase == 0 { + if ainvnm != 0 { + rcond = (1 / anorm) / ainvnm + } + return rcond + } + if kase == kase1 { + scale = impl.Dlatrs(uplo, blas.NoTrans, diag, normin, n, a, lda, work, work[2*n:]) + } else { + scale = impl.Dlatrs(uplo, blas.Trans, diag, normin, n, a, lda, work, work[2*n:]) + } + normin = true + if scale != 1 { + ix := bi.Idamax(n, work, 1) + xnorm := math.Abs(work[ix]) + if scale == 0 || scale < xnorm*smlnum { + return rcond + } + impl.Drscl(n, scale, work, 1) + } + } +} diff --git a/native/lapack_test.go b/native/lapack_test.go index 7c0f4fa4..cadc6ef5 100644 --- a/native/lapack_test.go +++ b/native/lapack_test.go @@ -104,6 +104,10 @@ func TestDrscl(t *testing.T) { testlapack.DrsclTest(t, impl) } +func TestDtrcon(t *testing.T) { + testlapack.DtrconTest(t, impl) +} + func TestIladlc(t *testing.T) { testlapack.IladlcTest(t, impl) } diff --git a/testlapack/dtrcon.go b/testlapack/dtrcon.go new file mode 100644 index 00000000..39ab34ee --- /dev/null +++ b/testlapack/dtrcon.go @@ -0,0 +1,179 @@ +package testlapack + +import ( + "math" + "math/rand" + "testing" + + "github.com/gonum/blas" + "github.com/gonum/floats" + "github.com/gonum/lapack" +) + +type Dtrconer interface { + Dgeconer + Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64 +} + +func DtrconTest(t *testing.T, impl Dtrconer) { + // Hand crafted tests. + for _, test := range []struct { + a []float64 + n int + uplo blas.Uplo + diag blas.Diag + condOne float64 + condInf float64 + }{ + { + a: []float64{ + 8, 5, 6, + 0, 7, 8, + 0, 0, 6, + }, + n: 3, + uplo: blas.Upper, + diag: blas.Unit, + condOne: 1.0 / 645, + condInf: 1.0 / 480, + }, + { + a: []float64{ + 8, 5, 6, + 0, 7, 8, + 0, 0, 6, + }, + n: 3, + uplo: blas.Upper, + diag: blas.NonUnit, + condOne: 0.137704918032787, + condInf: 0.157894736842105, + }, + { + a: []float64{ + 8, 0, 0, + 5, 7, 0, + 6, 8, 6, + }, + n: 3, + uplo: blas.Lower, + diag: blas.Unit, + condOne: 1.0 / 480, + condInf: 1.0 / 645, + }, + { + a: []float64{ + 8, 0, 0, + 5, 7, 0, + 6, 8, 6, + }, + n: 3, + uplo: blas.Lower, + diag: blas.NonUnit, + condOne: 0.157894736842105, + condInf: 0.137704918032787, + }, + } { + lda := test.n + work := make([]float64, 3*test.n) + for i := range work { + work[i] = rand.Float64() + } + iwork := make([]int, test.n) + for i := range iwork { + iwork[i] = rand.Int() + } + aCopy := make([]float64, len(test.a)) + copy(aCopy, test.a) + condOne := impl.Dtrcon(lapack.MaxColumnSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork) + if math.Abs(condOne-test.condOne) > 1e-14 { + t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne) + } + if !floats.Equal(aCopy, test.a) { + t.Errorf("a modified during call") + } + condInf := impl.Dtrcon(lapack.MaxRowSum, test.uplo, test.diag, test.n, test.a, lda, work, iwork) + if math.Abs(condInf-test.condInf) > 1e-14 { + t.Errorf("Inf norm mismatch. Want %v, got %v.", test.condInf, condInf) + } + if !floats.Equal(aCopy, test.a) { + t.Errorf("a modified during call") + } + } + + // Dtrcon does not match the Dgecon output in many cases. See + // https://github.com/xianyi/OpenBLAS/issues/636 + // TODO(btracey): Uncomment this when the mismatch between Dgecon and Dtrcon + // is understood. + /* + // Randomized tests against Dgecon. + for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} { + for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} { + for _, test := range []struct { + n, lda int + }{ + {3, 0}, + {4, 9}, + } { + for trial := 0; trial < 1; trial++ { + n := test.n + lda := test.lda + if lda == 0 { + lda = n + } + a := make([]float64, n*lda) + if trial == 0 { + for i := range a { + a[i] = float64(i + 2) + } + } else { + for i := range a { + a[i] = rand.NormFloat64() + } + } + + aDense := make([]float64, len(a)) + if uplo == blas.Upper { + for i := 0; i < n; i++ { + for j := i; j < n; j++ { + aDense[i*lda+j] = a[i*lda+j] + } + } + } else { + for i := 0; i < n; i++ { + for j := 0; j <= i; j++ { + aDense[i*lda+j] = a[i*lda+j] + } + } + } + if diag == blas.Unit { + for i := 0; i < n; i++ { + aDense[i*lda+i] = 1 + } + } + + ipiv := make([]int, n) + work := make([]float64, 4*n) + denseOne := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work) + denseInf := impl.Dlange(lapack.MaxRowSum, n, n, aDense, lda, work) + + aDenseLU := make([]float64, len(aDense)) + copy(aDenseLU, aDense) + impl.Dgetrf(n, n, aDenseLU, lda, ipiv) + iwork := make([]int, n) + want := impl.Dgecon(lapack.MaxColumnSum, n, aDenseLU, lda, denseOne, work, iwork) + got := impl.Dtrcon(lapack.MaxColumnSum, uplo, diag, n, a, lda, work, iwork) + if math.Abs(want-got) > 1e-14 { + t.Errorf("One norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got) + } + want = impl.Dgecon(lapack.MaxRowSum, n, aDenseLU, lda, denseInf, work, iwork) + got = impl.Dtrcon(lapack.MaxRowSum, uplo, diag, n, a, lda, work, iwork) + if math.Abs(want-got) > 1e-14 { + t.Errorf("Inf norm mismatch. Upper = %v, unit = %v, want %v, got %v", uplo == blas.Upper, diag == blas.Unit, want, got) + } + } + } + } + } + */ +}