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