blob: 8fba1848c660bf5c68c6f216aa08e8fa41c9096c [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" "math" "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/lapack" ) type Dpbconer interface { Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) float64 Dpbtrser } // DpbconTest tests Dpbcon by generating a random symmetric band matrix A and // checking that the estimated condition number is not too different from the // condition number computed via the explicit inverse of A. func DpbconTest(t *testing.T, impl Dpbconer) { rnd := rand.New(rand.NewSource(1)) for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} { 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 + 3} { dpbconTest(t, impl, uplo, n, kd, ldab, rnd) } } } } } func dpbconTest(t *testing.T, impl Dpbconer, uplo blas.Uplo, n, kd, ldab int, rnd *rand.Rand) { const ratioThresh = 10 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) } // Compute the norm of A. work := make([]float64, 3*n) aNorm := dlansb(lapack.MaxColumnSum, uplo, n, kd, ab, ldab, work) // Compute an estimate of rCond. iwork := make([]int, n) abFacCopy := make([]float64, len(abFac)) copy(abFacCopy, abFac) rCondGot := impl.Dpbcon(uplo, n, kd, abFac, ldab, aNorm, work, iwork) if !floats.Equal(abFac, abFacCopy) { t.Errorf("%v: unexpected modification of ab", name) } // Form the inverse of A to compute a good estimate of the condition number // rCondWant := 1/(norm(A) * norm(inv(A))) lda := max(1, n) aInv := make([]float64, n*lda) for i := 0; i < n; i++ { aInv[i*lda+i] = 1 } impl.Dpbtrs(uplo, n, kd, n, abFac, ldab, aInv, lda) aInvNorm := dlange(lapack.MaxColumnSum, n, n, aInv, lda) rCondWant := 1.0 if aNorm > 0 && aInvNorm > 0 { rCondWant = 1 / aNorm / aInvNorm } ratio := rCondTestRatio(rCondGot, rCondWant) if ratio >= ratioThresh { t.Errorf("%v: unexpected value of rcond. got=%v, want=%v (ratio=%v)", name, rCondGot, rCondWant, ratio) } } // rCondTestRatio returns a test ratio to compare two values of the reciprocal // of the condition number. // // This function corresponds to DGET06 in Reference LAPACK. func rCondTestRatio(rcond, rcondc float64) float64 { const eps = dlamchE switch { case rcond > 0 && rcondc > 0: return math.Max(rcond, rcondc)/math.Min(rcond, rcondc) - (1 - eps) case rcond > 0: return rcond / eps case rcondc > 0: return rcondc / eps default: return 0 } }