blob: f8a952aa6116ab1b857e32d025d80e91232c5158 [file] [log] [blame]
 // Copyright ©2015 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 ( "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" "gonum.org/v1/gonum/floats/scalar" ) type Dpotrfer interface { Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok bool) } func DpotrfTest(t *testing.T, impl Dpotrfer) { const tol = 1e-13 rnd := rand.New(rand.NewSource(1)) bi := blas64.Implementation() for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} { for tc, test := range []struct { n int lda int }{ {1, 0}, {2, 0}, {3, 0}, {10, 0}, {30, 0}, {63, 0}, {65, 0}, {127, 0}, {129, 0}, {500, 0}, {1, 10}, {2, 10}, {3, 10}, {10, 20}, {30, 50}, {63, 100}, {65, 100}, {127, 200}, {129, 200}, {500, 600}, } { n := test.n // Random diagonal matrix D with positive entries. d := make([]float64, n) Dlatm1(d, 4, 10000, false, 1, rnd) // Construct a positive definite matrix A as // A = U * D * Uᵀ // where U is a random orthogonal matrix. lda := test.lda if lda == 0 { lda = n } a := make([]float64, n*lda) Dlagsy(n, 0, d, a, lda, rnd, make([]float64, 2*n)) aCopy := make([]float64, len(a)) copy(aCopy, a) ok := impl.Dpotrf(uplo, n, a, lda) if !ok { t.Errorf("Case %v: unexpected failure for positive definite matrix", tc) continue } switch uplo { case blas.Upper: for i := 0; i < n; i++ { for j := 0; j < i; j++ { a[i*lda+j] = 0 } } case blas.Lower: for i := 0; i < n; i++ { for j := i + 1; j < n; j++ { a[i*lda+j] = 0 } } default: panic("bad uplo") } ans := make([]float64, len(a)) switch uplo { case blas.Upper: // Multiply Uᵀ * U. bi.Dsyrk(uplo, blas.Trans, n, n, 1, a, lda, 0, ans, lda) case blas.Lower: // Multiply L * Lᵀ. bi.Dsyrk(uplo, blas.NoTrans, n, n, 1, a, lda, 0, ans, lda) } match := true switch uplo { case blas.Upper: for i := 0; i < n; i++ { for j := i; j < n; j++ { if !scalar.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) { match = false } } } case blas.Lower: for i := 0; i < n; i++ { for j := 0; j <= i; j++ { if !scalar.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) { match = false } } } } if !match { t.Errorf("Case %v (uplo=%v,n=%v,lda=%v): unexpected result", tc, uplo, n, lda) } // Make one element of D negative so that A is not // positive definite, and check that Dpotrf fails. d[0] *= -1 Dlagsy(n, 0, d, a, lda, rnd, make([]float64, 2*n)) ok = impl.Dpotrf(uplo, n, a, lda) if ok { t.Errorf("Case %v: unexpected success for not positive definite matrix", tc) } } } }