blob: a08b6e728b08872068d3d276d00ca30aa9614396 [file] [log] [blame]
 // Copyright ©2019 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package testlapack import ( "fmt" "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" ) type Dpbtrfer interface { Dpbtrf(uplo blas.Uplo, n, kd int, ab []float64, ldab int) (ok bool) } // DpbtrfTest tests a band Cholesky factorization on random symmetric positive definite // band matrices by checking that the Cholesky factors multiply back to the original matrix. func DpbtrfTest(t *testing.T, impl Dpbtrfer) { // TODO(vladimir-ch): include expected-failure test case. // With the current implementation of Ilaenv the blocked code path is taken if kd > 64. // Unfortunately, with the block size nb=32 this also means that in Dpbtrf // it never happens that i2 <= 0 and the state coverage (unlike code coverage) is not complete. rnd := rand.New(rand.NewSource(1)) for _, n := range []int{0, 1, 2, 3, 4, 5, 64, 65, 66, 91, 96, 97, 101, 128, 130} { for _, kd := range []int{0, (n + 1) / 4, (3*n - 1) / 4, (5*n + 1) / 4} { for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} { for _, ldab := range []int{kd + 1, kd + 1 + 7} { dpbtrfTest(t, impl, uplo, n, kd, ldab, rnd) } } } } } func dpbtrfTest(t *testing.T, impl Dpbtrfer, uplo blas.Uplo, n, kd int, ldab int, rnd *rand.Rand) { const tol = 1e-12 name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,ldab=%v", string(uplo), n, kd, 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)) copy(abFac, ab) ok := impl.Dpbtrf(uplo, n, kd, abFac, ldab) if !ok { t.Fatalf("%v: bad test matrix, Dpbtrf failed", name) } // Reconstruct an symmetric band matrix from the Uᵀ*U or L*Lᵀ factorization, overwriting abFac. dsbmm(uplo, n, kd, abFac, ldab) // Compute and check the max-norm distance between the reconstructed and original matrix A. dist := distSymBand(uplo, n, kd, abFac, ldab, ab, ldab) if dist > tol*float64(n) { t.Errorf("%v: unexpected result, diff=%v", name, dist) } } // dsbmm computes a symmetric band matrix A // A = Uᵀ*U if uplo == blas.Upper, // A = L*Lᵀ if uplo == blas.Lower, // where U and L is an upper, respectively lower, triangular band matrix // stored on entry in ab. The result is stored in-place into ab. func dsbmm(uplo blas.Uplo, n, kd int, ab []float64, ldab int) { bi := blas64.Implementation() switch uplo { case blas.Upper: // Compute the product Uᵀ * U. for k := n - 1; k >= 0; k-- { klen := min(kd, n-k-1) // Number of stored off-diagonal elements in the row // Add a multiple of row k of the factor U to each of rows k+1 through n. if klen > 0 { bi.Dsyr(blas.Upper, klen, 1, ab[k*ldab+1:], 1, ab[(k+1)*ldab:], ldab-1) } // Scale row k by the diagonal element. bi.Dscal(klen+1, ab[k*ldab], ab[k*ldab:], 1) } case blas.Lower: // Compute the product L * Lᵀ. for k := n - 1; k >= 0; k-- { kc := max(0, kd-k) // Index of the first valid element in the row klen := kd - kc // Number of stored off-diagonal elements in the row // Compute the diagonal [k,k] element. ab[k*ldab+kd] = bi.Ddot(klen+1, ab[k*ldab+kc:], 1, ab[k*ldab+kc:], 1) // Compute the rest of column k. if klen > 0 { bi.Dtrmv(blas.Lower, blas.NoTrans, blas.NonUnit, klen, ab[(k-klen)*ldab+kd:], ldab-1, ab[k*ldab+kc:], 1) } } } }