lapack/testlapack: simplify randSymBand and use it in Dpb* tests

lapack/testlapack/dpbtf2.go

@@ -6,7 +6,6 @@ package testlapack
 
 import (
 	"fmt"
-	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
@@ -39,24 +38,8 @@ func dpbtf2Test(t *testing.T, impl Dpbtf2er, rnd *rand.Rand, uplo blas.Uplo, n,
 
 	name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,ldab=%v", string(uplo), n, kd, ldab)
 
-	// Allocate a band matrix and fill it with random numbers.
+	// Generate a random symmetric positive definite band matrix.
-	ab := make([]float64, n*ldab)
+	ab := randSymBand(uplo, n, kd, ldab, rnd)
-	for i := range ab {
-		ab[i] = rnd.NormFloat64()
-	}
-	// Make sure that the matrix U or L has a sufficiently positive diagonal.
-	switch uplo {
-	case blas.Upper:
-		for i := 0; i < n; i++ {
-			ab[i*ldab] = 2 + rnd.Float64()
-		}
-	case blas.Lower:
-		for i := 0; i < n; i++ {
-			ab[i*ldab+kd] = 2 + rnd.Float64()
-		}
-	}
-	// Compute U^T*U or L*L^T. The resulting (symmetric) matrix A will be positive definite.
-	dsbmm(uplo, n, kd, ab, ldab)
 
 	// Compute the Cholesky decomposition of A.
 	abFac := make([]float64, len(ab))
@@ -66,30 +49,12 @@ func dpbtf2Test(t *testing.T, impl Dpbtf2er, rnd *rand.Rand, uplo blas.Uplo, n,
 		t.Fatalf("%v: bad test matrix, Dpbtf2 failed", name)
 	}
 
-	if n == 0 {
-		return
-	}
-
 	// Reconstruct an symmetric band matrix from the U^T*U or L*L^T factorization, overwriting abFac.
 	dsbmm(uplo, n, kd, abFac, ldab)
 
 	// Compute and check the max-norm distance between the reconstructed and original matrix A.
-	var diff float64
+	dist := distSymBand(uplo, n, kd, abFac, ldab, ab, ldab)
-	switch uplo {
+	if dist > tol {
-	case blas.Upper:
+		t.Errorf("%v: unexpected result, diff=%v", name, dist)
-		for i := 0; i < n; i++ {
-			for j := 0; j < min(kd+1, n-i); j++ {
-				diff = math.Max(diff, math.Abs(abFac[i*ldab+j]-ab[i*ldab+j]))
-			}
-		}
-	case blas.Lower:
-		for i := 0; i < n; i++ {
-			for j := max(0, kd-i); j < kd+1; j++ {
-				diff = math.Max(diff, math.Abs(abFac[i*ldab+j]-ab[i*ldab+j]))
-			}
-		}
-	}
-	if diff > tol {
-		t.Errorf("%v: unexpected result, diff=%v", name, diff)
 	}
 }

lapack/testlapack/dpbtrf.go

@@ -6,7 +6,6 @@ package testlapack
 
 import (
 	"fmt"
-	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
@@ -44,24 +43,8 @@ func dpbtrfTest(t *testing.T, impl Dpbtrfer, uplo blas.Uplo, n, kd int, ldab int
 
 	name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,ldab=%v", string(uplo), n, kd, ldab)
 
-	// Allocate a band matrix and fill it with random numbers.
+	// Generate a random symmetric positive definite band matrix.
-	ab := make([]float64, n*ldab)
+	ab := randSymBand(uplo, n, kd, ldab, rnd)
-	for i := range ab {
-		ab[i] = rnd.NormFloat64()
-	}
-	// Make sure that the matrix U or L has a sufficiently positive diagonal.
-	switch uplo {
-	case blas.Upper:
-		for i := 0; i < n; i++ {
-			ab[i*ldab] = 2 + rnd.Float64()
-		}
-	case blas.Lower:
-		for i := 0; i < n; i++ {
-			ab[i*ldab+kd] = 2 + rnd.Float64()
-		}
-	}
-	// Compute U^T*U or L*L^T. The resulting (symmetric) matrix A will be positive definite.
-	dsbmm(uplo, n, kd, ab, ldab)
 
 	// Compute the Cholesky decomposition of A.
 	abFac := make([]float64, len(ab))
@@ -71,31 +54,13 @@ func dpbtrfTest(t *testing.T, impl Dpbtrfer, uplo blas.Uplo, n, kd int, ldab int
 		t.Fatalf("%v: bad test matrix, Dpbtrf failed", name)
 	}
 
-	if n == 0 {
-		return
-	}
-
 	// Reconstruct an symmetric band matrix from the U^T*U or L*L^T factorization, overwriting abFac.
 	dsbmm(uplo, n, kd, abFac, ldab)
 
 	// Compute and check the max-norm distance between the reconstructed and original matrix A.
-	var diff float64
+	dist := distSymBand(uplo, n, kd, abFac, ldab, ab, ldab)
-	switch uplo {
+	if dist > tol*float64(n) {
-	case blas.Upper:
+		t.Errorf("%v: unexpected result, diff=%v", name, dist)
-		for i := 0; i < n; i++ {
-			for j := 0; j < min(kd+1, n-i); j++ {
-				diff = math.Max(diff, math.Abs(abFac[i*ldab+j]-ab[i*ldab+j]))
-			}
-		}
-	case blas.Lower:
-		for i := 0; i < n; i++ {
-			for j := max(0, kd-i); j < kd+1; j++ {
-				diff = math.Max(diff, math.Abs(abFac[i*ldab+j]-ab[i*ldab+j]))
-			}
-		}
-	}
-	if diff > tol {
-		t.Errorf("%v: unexpected result, diff=%v", name, diff)
 	}
 }
 

lapack/testlapack/dpbtrs.go

@@ -46,25 +46,8 @@ func dpbtrsTest(t *testing.T, impl Dpbtrser, rnd *rand.Rand, uplo blas.Uplo, n,
 
 	name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,nrhs=%v,ldab=%v,ldb=%v", string(uplo), n, kd, nrhs, ldab, ldb)
 
-	// Allocate a band matrix and fill it with random numbers.
+	// Generate a random symmetric positive definite band matrix.
-	ab := make([]float64, n*ldab)
+	ab := randSymBand(uplo, n, kd, ldab, rnd)
-	for i := range ab {
-		ab[i] = rnd.NormFloat64()
-	}
-	// Make sure that the matrix U or L has a sufficiently positive diagonal.
-	switch uplo {
-	case blas.Upper:
-		for i := 0; i < n; i++ {
-			ab[i*ldab] = float64(n) + rnd.Float64()
-		}
-	case blas.Lower:
-		for i := 0; i < n; i++ {
-			ab[i*ldab+kd] = float64(n) + rnd.Float64()
-		}
-	}
-	// Compute U^T*U or L*L^T. The resulting (symmetric) matrix A will be
-	// positive definite and well-conditioned.
-	dsbmm(uplo, n, kd, ab, ldab)
 
 	// Compute the Cholesky decomposition of A.
 	abFac := make([]float64, len(ab))

lapack/testlapack/general.go

@@ -1031,57 +1031,49 @@ func equalApproxSymmetric(a, b blas64.Symmetric, tol float64) bool {
 }
 
 // randSymBand returns an n×n random symmetric positive definite band matrix
-// with kd diagonals, and the equivalent symmetric dense matrix.
+// with kd diagonals.
-func randSymBand(ul blas.Uplo, n, kd, ldab int, rnd *rand.Rand) (blas64.SymmetricBand, blas64.Symmetric) {
+func randSymBand(uplo blas.Uplo, n, kd, ldab int, rnd *rand.Rand) []float64 {
-	if n == 0 {
+	// Allocate a triangular band matrix U or L and fill it with random numbers.
-		return blas64.SymmetricBand{Uplo: ul, Stride: 1}, blas64.Symmetric{Uplo: ul, Stride: 1}
+	ab := make([]float64, n*ldab)
+	for i := range ab {
+		ab[i] = rnd.NormFloat64()
 	}
-	// A matrix is positive definite if and only if it has a Cholesky
+	// Make sure that the matrix U or L has a sufficiently positive diagonal.
-	// decomposition. Generate a random lower triangular band matrix L with
+	switch uplo {
-	// strictly positive diagonal, and construct the random symmetric band
+	case blas.Upper:
-	// matrix as L*L^T.
-	a := make([]float64, n*n)
 		for i := 0; i < n; i++ {
-		for j := max(0, i-kd); j <= i; j++ {
+			ab[i*ldab] = float64(n) + rnd.Float64()
-			a[i*n+j] = rnd.NormFloat64()
 		}
-		a[i*n+i] = math.Abs(a[i*n+i])
+	case blas.Lower:
-		// Add an extra amount to the diagonal in order to improve the condition number.
+		for i := 0; i < n; i++ {
-		a[i*n+i] += 1.5 * rnd.Float64()
+			ab[i*ldab+kd] = float64(n) + rnd.Float64()
 		}
-	agen := blas64.General{
-		Rows:   n,
-		Cols:   n,
-		Stride: n,
-		Data:   a,
 	}
+	// Compute U^T*U or L*L^T. The resulting (symmetric) matrix A will be
+	// positive definite and well-conditioned.
+	dsbmm(uplo, n, kd, ab, ldab)
+	return ab
+}
 
-	// Construct the SymDense from a*a^T
+// distSymBand returns the max-norm distance between the symmetric band matrices
-	c := make([]float64, n*n)
+// A and B.
-	cgen := blas64.General{
+func distSymBand(uplo blas.Uplo, n, kd int, a []float64, lda int, b []float64, ldb int) float64 {
-		Rows:   n,
+	var dist float64
-		Cols:   n,
+	switch uplo {
-		Stride: n,
+	case blas.Upper:
-		Data:   c,
+		for i := 0; i < n; i++ {
+			for j := 0; j < min(kd+1, n-i); j++ {
+				dist = math.Max(dist, math.Abs(a[i*lda+j]-b[i*ldb+j]))
 			}
-	blas64.Gemm(blas.NoTrans, blas.Trans, 1, agen, agen, 0, cgen)
-	sym := blas64.Symmetric{
-		N:      n,
-		Stride: n,
-		Data:   c,
-		Uplo:   ul,
 		}
-
+	case blas.Lower:
-	b := symToSymBand(ul, c, n, n, kd, ldab)
+		for i := 0; i < n; i++ {
-	band := blas64.SymmetricBand{
+			for j := max(0, kd-i); j < kd+1; j++ {
-		N:      n,
+				dist = math.Max(dist, math.Abs(a[i*lda+j]-b[i*ldb+j]))
-		K:      kd,
-		Stride: ldab,
-		Data:   b,
-		Uplo:   ul,
 			}
-
+		}
-	return band, sym
+	}
+	return dist
 }
 
 // symToSymBand takes the data in a Symmetric matrix and returns a