github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dpocon.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  	"log"
     9  	"testing"
    10  
    11  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  	"github.com/gopherd/gonum/floats/scalar"
    16  	"github.com/gopherd/gonum/lapack"
    17  )
    18  
    19  type Dpoconer interface {
    20  	Dpotrfer
    21  	Dgeconer
    22  	Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64
    23  	Dpocon(uplo blas.Uplo, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
    24  }
    25  
    26  func DpoconTest(t *testing.T, impl Dpoconer) {
    27  	for _, test := range []struct {
    28  		a    []float64
    29  		n    int
    30  		cond float64
    31  		uplo blas.Uplo
    32  	}{
    33  		{
    34  			a: []float64{
    35  				89, 59, 77,
    36  				0, 107, 59,
    37  				0, 0, 89,
    38  			},
    39  			uplo: blas.Upper,
    40  			n:    3,
    41  			cond: 0.050052137643379,
    42  		},
    43  		{
    44  			a: []float64{
    45  				89, 0, 0,
    46  				59, 107, 0,
    47  				77, 59, 89,
    48  			},
    49  			uplo: blas.Lower,
    50  			n:    3,
    51  			cond: 0.050052137643379,
    52  		},
    53  		// Dgecon does not match Dpocon for this case. https://github.com/xianyi/OpenBLAS/issues/664.
    54  		{
    55  			a: []float64{
    56  				2.9995576045549965, -2.0898894566158663, 3.965560740124006,
    57  				0, 1.9634729526261008, -2.8681002706874104,
    58  				0, 0, 5.502416670471008,
    59  			},
    60  			uplo: blas.Upper,
    61  			n:    3,
    62  			cond: 0.024054837369015203,
    63  		},
    64  	} {
    65  		n := test.n
    66  		a := make([]float64, len(test.a))
    67  		copy(a, test.a)
    68  		lda := n
    69  		uplo := test.uplo
    70  		work := make([]float64, 3*n)
    71  		anorm := impl.Dlansy(lapack.MaxColumnSum, uplo, n, a, lda, work)
    72  		// Compute cholesky decomposition
    73  		ok := impl.Dpotrf(uplo, n, a, lda)
    74  		if !ok {
    75  			t.Errorf("Bad test, matrix not positive definite")
    76  			continue
    77  		}
    78  		iwork := make([]int, n)
    79  		cond := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork)
    80  		// Error if not the same order, otherwise log the difference.
    81  		if !scalar.EqualWithinAbsOrRel(cond, test.cond, 1e0, 1e0) {
    82  			t.Errorf("Cond mismatch. Want %v, got %v.", test.cond, cond)
    83  		} else if !scalar.EqualWithinAbsOrRel(cond, test.cond, 1e-14, 1e-14) {
    84  			log.Printf("Dpocon cond mismatch. Want %v, got %v.", test.cond, cond)
    85  		}
    86  	}
    87  	rnd := rand.New(rand.NewSource(1))
    88  	bi := blas64.Implementation()
    89  	// Randomized tests compared against Dgecon.
    90  	for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
    91  		for _, test := range []struct {
    92  			n, lda int
    93  		}{
    94  			{3, 0},
    95  			{3, 5},
    96  		} {
    97  			for trial := 0; trial < 100; trial++ {
    98  				n := test.n
    99  				lda := test.lda
   100  				if lda == 0 {
   101  					lda = n
   102  				}
   103  				a := make([]float64, n*lda)
   104  				for i := range a {
   105  					a[i] = rnd.NormFloat64()
   106  				}
   107  
   108  				// Multiply a by itself to make it symmetric positive definite.
   109  				aCopy := make([]float64, len(a))
   110  				copy(aCopy, a)
   111  				bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, aCopy, lda, aCopy, lda, 0, a, lda)
   112  
   113  				aDat := make([]float64, len(aCopy))
   114  				copy(aDat, a)
   115  
   116  				aDense := make([]float64, len(a))
   117  				if uplo == blas.Upper {
   118  					for i := 0; i < n; i++ {
   119  						for j := i; j < n; j++ {
   120  							v := a[i*lda+j]
   121  							aDense[i*lda+j] = v
   122  							aDense[j*lda+i] = v
   123  						}
   124  					}
   125  				} else {
   126  					for i := 0; i < n; i++ {
   127  						for j := 0; j <= i; j++ {
   128  							v := a[i*lda+j]
   129  							aDense[i*lda+j] = v
   130  							aDense[j*lda+i] = v
   131  						}
   132  					}
   133  				}
   134  				work := make([]float64, 4*n)
   135  				iwork := make([]int, n)
   136  
   137  				anorm := impl.Dlansy(lapack.MaxColumnSum, uplo, n, a, lda, work)
   138  				ok := impl.Dpotrf(uplo, n, a, lda)
   139  				if !ok {
   140  					t.Errorf("Bad test, matrix not positive definite")
   141  					continue
   142  				}
   143  				got := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork)
   144  
   145  				denseNorm := impl.Dlange(lapack.MaxColumnSum, n, n, aDense, lda, work)
   146  				ipiv := make([]int, n)
   147  				impl.Dgetrf(n, n, aDense, lda, ipiv)
   148  				want := impl.Dgecon(lapack.MaxColumnSum, n, aDense, lda, denseNorm, work, iwork)
   149  				// Error if not the same order, otherwise log the difference.
   150  				if !scalar.EqualWithinAbsOrRel(want, got, 1e0, 1e0) {
   151  					t.Errorf("Dpocon and Dgecon mismatch. Dpocon %v, Dgecon %v.", got, want)
   152  				} else if !scalar.EqualWithinAbsOrRel(want, got, 1e-14, 1e-14) {
   153  					log.Printf("Dpocon and Dgecon mismatch. Dpocon %v, Dgecon %v.", got, want)
   154  				}
   155  			}
   156  		}
   157  	}
   158  }