gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/testlapack/dpttrf.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/lapack"
    15  )
    16  
    17  type Dpttrfer interface {
    18  	Dpttrf(n int, d, e []float64) (ok bool)
    19  }
    20  
    21  // DpttrfTest tests a tridiagonal Cholesky factorization on random symmetric
    22  // positive definite tridiagonal matrices by checking that the Cholesky factors
    23  // multiply back to the original matrix.
    24  func DpttrfTest(t *testing.T, impl Dpttrfer) {
    25  	rnd := rand.New(rand.NewSource(1))
    26  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50, 51, 52, 53, 54, 100} {
    27  		dpttrfTest(t, impl, rnd, n)
    28  	}
    29  }
    30  
    31  func dpttrfTest(t *testing.T, impl Dpttrfer, rnd *rand.Rand, n int) {
    32  	const tol = 1e-15
    33  
    34  	name := fmt.Sprintf("n=%v", n)
    35  
    36  	// Generate a random diagonally dominant symmetric tridiagonal matrix A.
    37  	d, e := newRandomSymTridiag(n, rnd)
    38  
    39  	// Make a copy of d and e to hold the factorization.
    40  	var dFac, eFac []float64
    41  	if n > 0 {
    42  		dFac = make([]float64, len(d))
    43  		copy(dFac, d)
    44  		if n > 1 {
    45  			eFac = make([]float64, len(e))
    46  			copy(eFac, e)
    47  		}
    48  	}
    49  
    50  	// Compute the Cholesky factorization of A.
    51  	ok := impl.Dpttrf(n, dFac, eFac)
    52  	if !ok {
    53  		t.Errorf("%v: bad test matrix, Dpttrf failed", name)
    54  		return
    55  	}
    56  
    57  	// Check the residual norm(L*D*Lᵀ - A)/(n * norm(A)).
    58  	resid := dpttrfResidual(n, d, e, dFac, eFac)
    59  	if resid > tol {
    60  		t.Errorf("%v: unexpected residual |L*D*Lᵀ - A|/(n * norm(A)); got %v, want <= %v", name, resid, tol)
    61  	}
    62  }
    63  
    64  func dpttrfResidual(n int, d, e, dFac, eFac []float64) float64 {
    65  	if n == 0 {
    66  		return 0
    67  	}
    68  
    69  	// Construct the difference L*D*Lᵀ - A.
    70  	dDiff := make([]float64, n)
    71  	eDiff := make([]float64, n-1)
    72  	dDiff[0] = dFac[0] - d[0]
    73  	for i, ef := range eFac {
    74  		de := dFac[i] * ef
    75  		dDiff[i+1] = de*ef + dFac[i+1] - d[i+1]
    76  		eDiff[i] = de - e[i]
    77  	}
    78  
    79  	// Compute the 1-norm of the difference L*D*Lᵀ - A.
    80  	var resid float64
    81  	if n == 1 {
    82  		resid = math.Abs(dDiff[0])
    83  	} else {
    84  		resid = math.Max(math.Abs(dDiff[0])+math.Abs(eDiff[0]), math.Abs(dDiff[n-1])+math.Abs(eDiff[n-2]))
    85  		for i := 1; i < n-1; i++ {
    86  			resid = math.Max(resid, math.Abs(dDiff[i])+math.Abs(eDiff[i-1])+math.Abs(eDiff[i]))
    87  		}
    88  	}
    89  
    90  	anorm := dlanst(lapack.MaxColumnSum, n, d, e)
    91  
    92  	// Compute norm(L*D*Lᵀ - A)/(n * norm(A)).
    93  	if anorm == 0 {
    94  		if resid != 0 {
    95  			return math.Inf(1)
    96  		}
    97  		return 0
    98  	}
    99  	return resid / float64(n) / anorm
   100  }
   101  
   102  func newRandomSymTridiag(n int, rnd *rand.Rand) (d, e []float64) {
   103  	if n == 0 {
   104  		return nil, nil
   105  	}
   106  
   107  	if n == 1 {
   108  		d = make([]float64, 1)
   109  		d[0] = rnd.Float64()
   110  		return d, nil
   111  	}
   112  
   113  	// Allocate the diagonal d and fill it with numbers from [0,1).
   114  	d = make([]float64, n)
   115  	dlarnv(d, 1, rnd)
   116  	// Allocate the subdiagonal e and fill it with numbers from [-1,1).
   117  	e = make([]float64, n-1)
   118  	dlarnv(e, 2, rnd)
   119  
   120  	// Make A diagonally dominant by adding the absolute value of off-diagonal
   121  	// elements to it.
   122  	d[0] += math.Abs(e[0])
   123  	for i := 1; i < n-1; i++ {
   124  		d[i] += math.Abs(e[i]) + math.Abs(e[i-1])
   125  	}
   126  	d[n-1] += math.Abs(e[n-2])
   127  
   128  	return d, e
   129  }