lapack/testlapack: rewrite test for Dgecon from table-based to randomized

lapack/testlapack/dgecon.go

@@ -5,92 +5,87 @@
 package testlapack
 
 import (
+	"fmt"
 	"testing"
 
+	"golang.org/x/exp/rand"
+
 	"gonum.org/v1/gonum/floats"
 	"gonum.org/v1/gonum/lapack"
 )
 
 type Dgeconer interface {
-	Dlanger
-	Dgetrfer
 	Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
+
+	Dgetrier
+	Dlanger
 }
 
 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,
-		},
-		// Dgecon does not match Dpocon for this case. https://github.com/xianyi/OpenBLAS/issues/664.
-		{
-			a: []float64{
-				2.9995576045549965, -2.0898894566158663, 3.965560740124006,
-				-2.0898894566158663, 1.9634729526261008, -2.8681002706874104,
-				3.965560740124006, -2.8681002706874104, 5.502416670471008,
-			},
-			m:       3,
-			n:       3,
-			condOne: 0.024054837369015203,
-			condInf: 0.024054837369015203,
-		},
-	} {
-		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)
-
-		// Error if not the same order, otherwise log the difference.
-		if !floats.EqualWithinAbsOrRel(condOne, test.condOne, 1e0, 1e0) {
-			t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne)
-		} else if !floats.EqualWithinAbsOrRel(condOne, test.condOne, 1e-14, 1e-14) {
-			t.Logf("Dgecon one norm mismatch. Want %v, got %v.", test.condOne, condOne)
-		}
-		if !floats.EqualWithinAbsOrRel(condInf, test.condInf, 1e0, 1e0) {
-			t.Errorf("One norm mismatch. Want %v, got %v.", test.condInf, condInf)
-		} else if !floats.EqualWithinAbsOrRel(condInf, test.condInf, 1e-14, 1e-14) {
-			t.Logf("Dgecon one norm mismatch. Want %v, got %v.", test.condInf, condInf)
+	rnd := rand.New(rand.NewSource(1))
+	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
+		for _, lda := range []int{max(1, n), n + 3} {
+			dgeconTest(t, impl, rnd, n, lda)
+		}
+	}
+}
+
+func dgeconTest(t *testing.T, impl Dgeconer, rnd *rand.Rand, n, lda int) {
+	const ratioThresh = 10
+
+	// Generate a random square matrix A with elements uniformly in [-1,1).
+	a := make([]float64, max(0, (n-1)*lda+n))
+	for i := range a {
+		a[i] = 2*rnd.Float64() - 1
+	}
+
+	// Allocate work slices.
+	iwork := make([]int, n)
+	work := make([]float64, max(1, 4*n))
+
+	// Compute the LU factorization of A.
+	aFac := make([]float64, len(a))
+	copy(aFac, a)
+	ipiv := make([]int, n)
+	ok := impl.Dgetrf(n, n, aFac, lda, ipiv)
+	if !ok {
+		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetrf failed", n, lda)
+	}
+	aFacCopy := make([]float64, len(aFac))
+	copy(aFacCopy, aFac)
+
+	// Compute the inverse A^{-1} from the LU factorization.
+	aInv := make([]float64, len(aFac))
+	copy(aInv, aFac)
+	ok = impl.Dgetri(n, aInv, lda, ipiv, work, len(work))
+	if !ok {
+		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetri failed", n, lda)
+	}
+
+	for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
+		name := fmt.Sprintf("norm%v,n=%v,lda=%v", string(norm), n, lda)
+
+		// Compute the norm of A and A^{-1}.
+		aNorm := impl.Dlange(norm, n, n, a, lda, work)
+		aInvNorm := impl.Dlange(norm, n, n, aInv, lda, work)
+
+		// Compute a good estimate of the condition number
+		//  rcondWant := 1/(norm(A) * norm(inv(A)))
+		rcondWant := 1.0
+		if aNorm > 0 && aInvNorm > 0 {
+			rcondWant = 1 / aNorm / aInvNorm
+		}
+
+		// Compute an estimate of rcond using the LU factorization and Dgecon.
+		rcondGot := impl.Dgecon(norm, n, aFac, lda, aNorm, work, iwork)
+		if !floats.Equal(aFac, aFacCopy) {
+			t.Errorf("%v: unexpected modification of aFac", name)
+		}
+
+		ratio := rCondTestRatio(rcondGot, rcondWant)
+		if ratio >= ratioThresh {
+			t.Errorf("%v: unexpected value of rcond; got=%v, want=%v (ratio=%v)",
+				name, rcondGot, rcondWant, ratio)
 		}
 	}
 }