github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dpbcon.go (about)

     1  // Copyright ©2019 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  	"fmt"
     9  	"math"
    10  	"testing"
    11  
    12  	"math/rand"
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/floats"
    15  	"github.com/gopherd/gonum/lapack"
    16  )
    17  
    18  type Dpbconer interface {
    19  	Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) float64
    20  
    21  	Dpbtrser
    22  }
    23  
    24  // DpbconTest tests Dpbcon by generating a random symmetric band matrix A and
    25  // checking that the estimated condition number is not too different from the
    26  // condition number computed via the explicit inverse of A.
    27  func DpbconTest(t *testing.T, impl Dpbconer) {
    28  	rnd := rand.New(rand.NewSource(1))
    29  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
    30  		for _, kd := range []int{0, (n + 1) / 4, (3*n - 1) / 4, (5*n + 1) / 4} {
    31  			for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    32  				for _, ldab := range []int{kd + 1, kd + 1 + 3} {
    33  					dpbconTest(t, impl, uplo, n, kd, ldab, rnd)
    34  				}
    35  			}
    36  		}
    37  	}
    38  }
    39  
    40  func dpbconTest(t *testing.T, impl Dpbconer, uplo blas.Uplo, n, kd, ldab int, rnd *rand.Rand) {
    41  	const ratioThresh = 10
    42  
    43  	name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,ldab=%v", string(uplo), n, kd, ldab)
    44  
    45  	// Generate a random symmetric positive definite band matrix.
    46  	ab := randSymBand(uplo, n, kd, ldab, rnd)
    47  
    48  	// Compute the Cholesky decomposition of A.
    49  	abFac := make([]float64, len(ab))
    50  	copy(abFac, ab)
    51  	ok := impl.Dpbtrf(uplo, n, kd, abFac, ldab)
    52  	if !ok {
    53  		t.Fatalf("%v: bad test matrix, Dpbtrf failed", name)
    54  	}
    55  
    56  	// Compute the norm of A.
    57  	work := make([]float64, 3*n)
    58  	aNorm := dlansb(lapack.MaxColumnSum, uplo, n, kd, ab, ldab, work)
    59  
    60  	// Compute an estimate of rCond.
    61  	iwork := make([]int, n)
    62  	abFacCopy := make([]float64, len(abFac))
    63  	copy(abFacCopy, abFac)
    64  	rCondGot := impl.Dpbcon(uplo, n, kd, abFac, ldab, aNorm, work, iwork)
    65  
    66  	if !floats.Equal(abFac, abFacCopy) {
    67  		t.Errorf("%v: unexpected modification of ab", name)
    68  	}
    69  
    70  	// Form the inverse of A to compute a good estimate of the condition number
    71  	//  rCondWant := 1/(norm(A) * norm(inv(A)))
    72  	lda := max(1, n)
    73  	aInv := make([]float64, n*lda)
    74  	for i := 0; i < n; i++ {
    75  		aInv[i*lda+i] = 1
    76  	}
    77  	impl.Dpbtrs(uplo, n, kd, n, abFac, ldab, aInv, lda)
    78  	aInvNorm := dlange(lapack.MaxColumnSum, n, n, aInv, lda)
    79  	rCondWant := 1.0
    80  	if aNorm > 0 && aInvNorm > 0 {
    81  		rCondWant = 1 / aNorm / aInvNorm
    82  	}
    83  
    84  	ratio := rCondTestRatio(rCondGot, rCondWant)
    85  	if ratio >= ratioThresh {
    86  		t.Errorf("%v: unexpected value of rcond. got=%v, want=%v (ratio=%v)", name, rCondGot, rCondWant, ratio)
    87  	}
    88  }
    89  
    90  // rCondTestRatio returns a test ratio to compare two values of the reciprocal
    91  // of the condition number.
    92  //
    93  // This function corresponds to DGET06 in Reference LAPACK.
    94  func rCondTestRatio(rcond, rcondc float64) float64 {
    95  	const eps = dlamchE
    96  	switch {
    97  	case rcond > 0 && rcondc > 0:
    98  		return math.Max(rcond, rcondc)/math.Min(rcond, rcondc) - (1 - eps)
    99  	case rcond > 0:
   100  		return rcond / eps
   101  	case rcondc > 0:
   102  		return rcondc / eps
   103  	default:
   104  		return 0
   105  	}
   106  }