Add Dpocon and tests

cgo/lapack.go

@@ -524,6 +524,31 @@ 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)
+	}
+	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.
 // The condition number computed may be based on the 1-norm or the ∞-norm.
 //

cgo/lapack_test.go

@@ -87,6 +87,13 @@ 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)
 }

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)
+		}
+	}
+}

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)
 }

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)
+				}
+			}
+		}
+	}
+}