gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/testlapack/dptcon.go (about)

     1  // Copyright ©2023 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  	"golang.org/x/exp/rand"
    13  
    14  	"gonum.org/v1/gonum/floats"
    15  	"gonum.org/v1/gonum/lapack"
    16  )
    17  
    18  type Dptconer interface {
    19  	Dptcon(n int, d, e []float64, anorm float64, work []float64) (rcond float64)
    20  
    21  	Dpttrf(n int, d, e []float64) (ok bool)
    22  	Dpttrs(n, nrhs int, d, e []float64, b []float64, ldb int)
    23  }
    24  
    25  func DptconTest(t *testing.T, impl Dptconer) {
    26  	rnd := rand.New(rand.NewSource(1))
    27  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50, 51, 52, 53, 54, 100} {
    28  		dptconTest(t, impl, rnd, n)
    29  	}
    30  }
    31  
    32  func dptconTest(t *testing.T, impl Dptconer, rnd *rand.Rand, n int) {
    33  	const tol = 1e-15
    34  
    35  	name := fmt.Sprintf("n=%v", n)
    36  
    37  	// Generate a random diagonally dominant symmetric tridiagonal matrix A.
    38  	d, e := newRandomSymTridiag(n, rnd)
    39  	aNorm := dlanst(lapack.MaxColumnSum, n, d, e)
    40  
    41  	// Compute the Cholesky factorization of A.
    42  	ok := impl.Dpttrf(n, d, e)
    43  	if !ok {
    44  		t.Errorf("%v: bad test matrix, Dpttrf failed", name)
    45  		return
    46  	}
    47  
    48  	// Compute the reciprocal of the condition number of A.
    49  	dCopy := make([]float64, len(d))
    50  	copy(dCopy, d)
    51  	eCopy := make([]float64, len(e))
    52  	copy(eCopy, e)
    53  	work := make([]float64, 3*n)
    54  	rcondGot := impl.Dptcon(n, d, e, aNorm, work)
    55  
    56  	// Check that Dptcon didn't modify d and e.
    57  	if !floats.Equal(d, dCopy) {
    58  		t.Errorf("%v: unexpected modification of d", name)
    59  	}
    60  	if !floats.Equal(e, eCopy) {
    61  		t.Errorf("%v: unexpected modification of e", name)
    62  	}
    63  
    64  	// Compute the norm of A⁻¹.
    65  	aInv, lda := make([]float64, n*n), max(1, n)
    66  	for i := 0; i < n; i++ {
    67  		aInv[i*lda+i] = 1
    68  	}
    69  	impl.Dpttrs(n, n, d, e, aInv, lda)
    70  	aInvNorm := dlange(lapack.MaxColumnSum, n, n, aInv, lda)
    71  
    72  	rcondWant := 1.0
    73  	if aNorm > 0 && aInvNorm > 0 {
    74  		rcondWant = 1 / aNorm / aInvNorm
    75  	}
    76  
    77  	diff := math.Abs(rcondGot - rcondWant)
    78  	if diff > tol {
    79  		t.Errorf("%v: unexpected value of rcond. got=%v, want=%v (diff=%v)", name, rcondGot, rcondWant, diff)
    80  	}
    81  }