From ce6986a6783167c78406a5bbd525cedef6413c1c Mon Sep 17 00:00:00 2001 From: Vladimir Chalupecky Date: Sun, 16 Jun 2019 18:51:40 +0200 Subject: [PATCH] lapack/testlapack: simplify randSymBand and use it in Dpb* tests --- lapack/testlapack/dpbtf2.go | 45 +++---------------- lapack/testlapack/dpbtrf.go | 45 +++---------------- lapack/testlapack/dpbtrs.go | 21 +-------- lapack/testlapack/general.go | 84 ++++++++++++++++-------------------- 4 files changed, 50 insertions(+), 145 deletions(-) diff --git a/lapack/testlapack/dpbtf2.go b/lapack/testlapack/dpbtf2.go index c882119c..3c317225 100644 --- a/lapack/testlapack/dpbtf2.go +++ b/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. - ab := make([]float64, n*ldab) - 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) + // Generate a random symmetric positive definite band matrix. + ab := randSymBand(uplo, n, kd, ldab, rnd) // 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 - switch uplo { - case blas.Upper: - 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) + dist := distSymBand(uplo, n, kd, abFac, ldab, ab, ldab) + if dist > tol { + t.Errorf("%v: unexpected result, diff=%v", name, dist) } } diff --git a/lapack/testlapack/dpbtrf.go b/lapack/testlapack/dpbtrf.go index 89ec30fb..0125735e 100644 --- a/lapack/testlapack/dpbtrf.go +++ b/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. - ab := make([]float64, n*ldab) - 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) + // Generate a random symmetric positive definite band matrix. + ab := randSymBand(uplo, n, kd, ldab, rnd) // 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 - switch uplo { - case blas.Upper: - 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) + dist := distSymBand(uplo, n, kd, abFac, ldab, ab, ldab) + if dist > tol*float64(n) { + t.Errorf("%v: unexpected result, diff=%v", name, dist) } } diff --git a/lapack/testlapack/dpbtrs.go b/lapack/testlapack/dpbtrs.go index b68cbd65..b9daa86e 100644 --- a/lapack/testlapack/dpbtrs.go +++ b/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. - ab := make([]float64, n*ldab) - 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) + // Generate a random symmetric positive definite band matrix. + ab := randSymBand(uplo, n, kd, ldab, rnd) // Compute the Cholesky decomposition of A. abFac := make([]float64, len(ab)) diff --git a/lapack/testlapack/general.go b/lapack/testlapack/general.go index da23eb8a..07116637 100644 --- a/lapack/testlapack/general.go +++ b/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. -func randSymBand(ul blas.Uplo, n, kd, ldab int, rnd *rand.Rand) (blas64.SymmetricBand, blas64.Symmetric) { - if n == 0 { - return blas64.SymmetricBand{Uplo: ul, Stride: 1}, blas64.Symmetric{Uplo: ul, Stride: 1} +// with kd diagonals. +func randSymBand(uplo blas.Uplo, n, kd, ldab int, rnd *rand.Rand) []float64 { + // Allocate a triangular band matrix U or L and fill it with random numbers. + 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 - // decomposition. Generate a random lower triangular band matrix L with - // strictly positive diagonal, and construct the random symmetric band - // 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++ { - a[i*n+j] = 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() } - a[i*n+i] = math.Abs(a[i*n+i]) - // Add an extra amount to the diagonal in order to improve the condition number. - a[i*n+i] += 1.5 * 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 - c := make([]float64, n*n) - cgen := blas64.General{ - Rows: n, - Cols: n, - Stride: n, - Data: c, +// distSymBand returns the max-norm distance between the symmetric band matrices +// A and B. +func distSymBand(uplo blas.Uplo, n, kd int, a []float64, lda int, b []float64, ldb int) float64 { + var dist float64 + switch uplo { + case blas.Upper: + 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])) + } + } + case blas.Lower: + for i := 0; i < n; i++ { + for j := max(0, kd-i); j < kd+1; 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, - } - - b := symToSymBand(ul, c, n, n, kd, ldab) - band := blas64.SymmetricBand{ - N: n, - K: kd, - Stride: ldab, - Data: b, - Uplo: ul, - } - - return band, sym + return dist } // symToSymBand takes the data in a Symmetric matrix and returns a