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