Add Dgecon and tests, plus small fix to Dlacn2

cgo/lapack.go

@@ -129,6 +129,33 @@ func (impl Implementation) Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok
 	return clapack.Dpotrf(ul, n, a, lda)
 }
 
+// Dgecon estimates the reciprocal of the condition number of the n×n matrix A
+// given the LU decomposition of the matrix. The condition number computed may
+// be based on the 1-norm or the ∞-norm.
+//
+// The slice a contains the result of the LU decomposition of A as computed by Dgetrf.
+//
+// anorm is the corresponding 1-norm or ∞-norm of the original matrix A.
+//
+// work is a temporary data slice of length at least 4*n and Dgecon will panic otherwise.
+//
+// iwork is a temporary data slice of length at least n and Dgecon will panic otherwise.
+func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 {
+	checkMatrix(n, n, a, lda)
+	if norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum {
+		panic("bad norm")
+	}
+	if len(work) < 4*n {
+		panic(badWork)
+	}
+	if len(iwork) < n {
+		panic(badWork)
+	}
+	rcond := make([]float64, 1)
+	clapack.Dgecon(byte(norm), n, a, lda, anorm, rcond)
+	return rcond[0]
+}
+
 // Dgelq2 computes the LQ factorization of the m×n matrix A.
 //
 // In an LQ factorization, L is a lower triangular m×n matrix, and Q is an n×n

cgo/lapack_test.go

@@ -28,6 +28,10 @@ func TestDpotrf(t *testing.T) {
 	testlapack.DpotrfTest(t, impl)
 }
 
+func TestDgecon(t *testing.T) {
+	testlapack.DgeconTest(t, impl)
+}
+
 func TestDgelq2(t *testing.T) {
 	testlapack.Dgelq2Test(t, impl)
 }

native/dgecon.go

@@ -0,0 +1,77 @@
+package native
+
+import (
+	"math"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+	"github.com/gonum/lapack"
+)
+
+// Dgecon estimates the reciprocal of the condition number of the n×n matrix A
+// given the LU decomposition of the matrix. The condition number computed may
+// be based on the 1-norm or the ∞-norm.
+//
+// The slice a contains the result of the LU decomposition of A as computed by Dgetrf.
+//
+// anorm is the corresponding 1-norm or ∞-norm of the original matrix A.
+//
+// work is a temporary data slice of length at least 4*n and Dgecon will panic otherwise.
+//
+// iwork is a temporary data slice of length at least n and Dgecon will panic otherwise.
+func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 {
+	checkMatrix(n, n, a, lda)
+	if norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum {
+		panic(badNorm)
+	}
+	if len(work) < 4*n {
+		panic(badWork)
+	}
+	if len(iwork) < n {
+		panic(badWork)
+	}
+
+	if n == 0 {
+		return 0
+	} else if anorm == 0 {
+		return 1
+	}
+
+	bi := blas64.Implementation()
+	var rcond, ainvnm float64
+	var kase int
+	var normin bool
+	isave := new([3]int)
+	onenrm := norm == lapack.MaxColumnSum
+	smlnum := dlamchS
+	kase1 := 2
+	if onenrm {
+		kase1 = 1
+	}
+	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
+		}
+		var sl, su float64
+		if kase == kase1 {
+			sl = impl.Dlatrs(blas.Lower, blas.NoTrans, blas.Unit, normin, n, a, lda, work, work[2*n:])
+			su = impl.Dlatrs(blas.Upper, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[3*n:])
+		} else {
+			su = impl.Dlatrs(blas.Upper, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[3*n:])
+			sl = impl.Dlatrs(blas.Lower, blas.Trans, blas.Unit, normin, n, a, lda, work, work[2*n:])
+		}
+		scale := sl * su
+		normin = true
+		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/dlacn2.go

@@ -24,7 +24,10 @@ func (impl Implementation) Dlacn2(n int, v, x []float64, isgn []int, est float64
 	if len(isgn) < n {
 		panic("lapack: insufficient isgn length")
 	}
-	if isave[0] < 1 || isave[0] > 5 {
+	if isave[0] < 0 || isave[0] > 5 {
+		panic("lapack: bad isave value")
+	}
+	if isave[0] == 0 && kase != 0 {
 		panic("lapack: bad isave value")
 	}
 	itmax := 5
@@ -39,7 +42,7 @@ func (impl Implementation) Dlacn2(n int, v, x []float64, isgn []int, est float64
 	}
 	switch isave[0] {
 	default:
-		panic("unknown case")
+		panic("unreachable")
 	case 1:
 		if n == 1 {
 			v[0] = x[0]

native/lapack_test.go

@@ -12,6 +12,10 @@ import (
 
 var impl = Implementation{}
 
+func TestDgecon(t *testing.T) {
+	testlapack.DgeconTest(t, impl)
+}
+
 func TestDgelqf(t *testing.T) {
 	testlapack.DgelqfTest(t, impl)
 }

testlapack/dgecon.go

@@ -0,0 +1,74 @@
+package testlapack
+
+import (
+	"math"
+	"testing"
+
+	"github.com/gonum/lapack"
+)
+
+type Dgeconer interface {
+	Dlanger
+	Dgetrfer
+	Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
+}
+
+func DgeconTest(t *testing.T, impl Dgeconer) {
+	for _, test := range []struct {
+		m       int
+		n       int
+		a       []float64
+		condOne float64
+		condInf float64
+	}{
+		{
+			a: []float64{
+				8, 1, 6,
+				3, 5, 7,
+				4, 9, 2,
+			},
+			m:       3,
+			n:       3,
+			condOne: 3.0 / 16,
+			condInf: 3.0 / 16,
+		},
+		{
+			a: []float64{
+				2, 9, 3, 2,
+				10, 9, 9, 3,
+				1, 1, 5, 2,
+				8, 4, 10, 2,
+			},
+			m:       4,
+			n:       4,
+			condOne: 0.024740155174938,
+			condInf: 0.012034465570035,
+		},
+	} {
+		m := test.m
+		n := test.n
+		lda := n
+		a := make([]float64, len(test.a))
+		copy(a, test.a)
+		ipiv := make([]int, min(m, n))
+
+		// Find the norms of the original matrix.
+		work := make([]float64, 4*n)
+		oneNorm := impl.Dlange(lapack.MaxColumnSum, m, n, a, lda, work)
+		infNorm := impl.Dlange(lapack.MaxRowSum, m, n, a, lda, work)
+
+		// Compute LU factorization of a.
+		impl.Dgetrf(m, n, a, lda, ipiv)
+
+		// Compute the condition number
+		iwork := make([]int, n)
+		condOne := impl.Dgecon(lapack.MaxColumnSum, n, a, lda, oneNorm, work, iwork)
+		condInf := impl.Dgecon(lapack.MaxRowSum, n, a, lda, infNorm, work, iwork)
+		if math.Abs(condOne-test.condOne) > 1e-13 {
+			t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne)
+		}
+		if math.Abs(condInf-test.condInf) > 1e-13 {
+			t.Errorf("Inf norm mismatch. Want %v, got %v.", test.condInf, condInf)
+		}
+	}
+}