gonum.org/v1/gonum@v0.14.0/lapack/testlapack/dpotrf.go (about)

     1  // Copyright ©2015 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package testlapack
     6  
     7  import (
     8  	"testing"
     9  
    10  	"golang.org/x/exp/rand"
    11  
    12  	"gonum.org/v1/gonum/blas"
    13  	"gonum.org/v1/gonum/blas/blas64"
    14  	"gonum.org/v1/gonum/floats/scalar"
    15  )
    16  
    17  type Dpotrfer interface {
    18  	Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok bool)
    19  }
    20  
    21  func DpotrfTest(t *testing.T, impl Dpotrfer) {
    22  	const tol = 1e-13
    23  	rnd := rand.New(rand.NewSource(1))
    24  	bi := blas64.Implementation()
    25  	for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    26  		for tc, test := range []struct {
    27  			n   int
    28  			lda int
    29  		}{
    30  			{1, 0},
    31  			{2, 0},
    32  			{3, 0},
    33  			{10, 0},
    34  			{30, 0},
    35  			{63, 0},
    36  			{65, 0},
    37  			{127, 0},
    38  			{129, 0},
    39  			{500, 0},
    40  			{1, 10},
    41  			{2, 10},
    42  			{3, 10},
    43  			{10, 20},
    44  			{30, 50},
    45  			{63, 100},
    46  			{65, 100},
    47  			{127, 200},
    48  			{129, 200},
    49  			{500, 600},
    50  		} {
    51  			n := test.n
    52  
    53  			// Random diagonal matrix D with positive entries.
    54  			d := make([]float64, n)
    55  			Dlatm1(d, 4, 10000, false, 1, rnd)
    56  
    57  			// Construct a positive definite matrix A as
    58  			//  A = U * D * Uᵀ
    59  			// where U is a random orthogonal matrix.
    60  			lda := test.lda
    61  			if lda == 0 {
    62  				lda = n
    63  			}
    64  			a := make([]float64, n*lda)
    65  			Dlagsy(n, 0, d, a, lda, rnd, make([]float64, 2*n))
    66  
    67  			aCopy := make([]float64, len(a))
    68  			copy(aCopy, a)
    69  
    70  			ok := impl.Dpotrf(uplo, n, a, lda)
    71  			if !ok {
    72  				t.Errorf("Case %v: unexpected failure for positive definite matrix", tc)
    73  				continue
    74  			}
    75  
    76  			switch uplo {
    77  			case blas.Upper:
    78  				for i := 0; i < n; i++ {
    79  					for j := 0; j < i; j++ {
    80  						a[i*lda+j] = 0
    81  					}
    82  				}
    83  			case blas.Lower:
    84  				for i := 0; i < n; i++ {
    85  					for j := i + 1; j < n; j++ {
    86  						a[i*lda+j] = 0
    87  					}
    88  				}
    89  			default:
    90  				panic("bad uplo")
    91  			}
    92  
    93  			ans := make([]float64, len(a))
    94  			switch uplo {
    95  			case blas.Upper:
    96  				// Multiply Uᵀ * U.
    97  				bi.Dsyrk(uplo, blas.Trans, n, n, 1, a, lda, 0, ans, lda)
    98  			case blas.Lower:
    99  				// Multiply L * Lᵀ.
   100  				bi.Dsyrk(uplo, blas.NoTrans, n, n, 1, a, lda, 0, ans, lda)
   101  			}
   102  
   103  			match := true
   104  			switch uplo {
   105  			case blas.Upper:
   106  				for i := 0; i < n; i++ {
   107  					for j := i; j < n; j++ {
   108  						if !scalar.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
   109  							match = false
   110  						}
   111  					}
   112  				}
   113  			case blas.Lower:
   114  				for i := 0; i < n; i++ {
   115  					for j := 0; j <= i; j++ {
   116  						if !scalar.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
   117  							match = false
   118  						}
   119  					}
   120  				}
   121  			}
   122  			if !match {
   123  				t.Errorf("Case %v (uplo=%v,n=%v,lda=%v): unexpected result", tc, uplo, n, lda)
   124  			}
   125  
   126  			// Make one element of D negative so that A is not
   127  			// positive definite, and check that Dpotrf fails.
   128  			d[0] *= -1
   129  			Dlagsy(n, 0, d, a, lda, rnd, make([]float64, 2*n))
   130  			ok = impl.Dpotrf(uplo, n, a, lda)
   131  			if ok {
   132  				t.Errorf("Case %v: unexpected success for not positive definite matrix", tc)
   133  			}
   134  		}
   135  	}
   136  }