gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/testlapack/dpttrs.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  	"testing"
    10  
    11  	"golang.org/x/exp/rand"
    12  	"gonum.org/v1/gonum/blas/blas64"
    13  	"gonum.org/v1/gonum/lapack"
    14  )
    15  
    16  type Dpttrser interface {
    17  	Dpttrs(n, nrhs int, d, e []float64, b []float64, ldb int)
    18  
    19  	Dpttrfer
    20  }
    21  
    22  func DpttrsTest(t *testing.T, impl Dpttrser) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50, 51, 52, 53, 54, 100} {
    25  		for _, nrhs := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50} {
    26  			for _, ldb := range []int{max(1, nrhs), nrhs + 3} {
    27  				dpttrsTest(t, impl, rnd, n, nrhs, ldb)
    28  			}
    29  		}
    30  	}
    31  }
    32  
    33  func dpttrsTest(t *testing.T, impl Dpttrser, rnd *rand.Rand, n, nrhs, ldb int) {
    34  	const tol = 1e-15
    35  
    36  	name := fmt.Sprintf("n=%v", n)
    37  
    38  	// Generate a random diagonally dominant symmetric tridiagonal matrix A.
    39  	d, e := newRandomSymTridiag(n, rnd)
    40  
    41  	// Make a copy of d and e to hold the factorization.
    42  	dFac := make([]float64, len(d))
    43  	copy(dFac, d)
    44  	eFac := make([]float64, len(e))
    45  	copy(eFac, e)
    46  
    47  	// Compute the Cholesky factorization of A.
    48  	ok := impl.Dpttrf(n, dFac, eFac)
    49  	if !ok {
    50  		t.Errorf("%v: bad test matrix, Dpttrf failed", name)
    51  		return
    52  	}
    53  
    54  	// Generate a random solution matrix X.
    55  	xWant := randomGeneral(n, nrhs, ldb, rnd)
    56  
    57  	// Compute the right-hand side.
    58  	b := zeros(n, nrhs, ldb)
    59  	dstmm(n, nrhs, d, e, xWant.Data, xWant.Stride, b.Data, b.Stride)
    60  
    61  	// Solve A*X=B.
    62  	impl.Dpttrs(n, nrhs, dFac, eFac, b.Data, b.Stride)
    63  
    64  	resid := dpttrsResidual(b, xWant)
    65  	if resid > tol {
    66  		t.Errorf("%v: unexpected solution: |diff| = %v, want <= %v", name, resid, tol)
    67  	}
    68  }
    69  
    70  // dstmm computes the matrix-matrix product
    71  //
    72  //	C = A*B
    73  //
    74  // where A is an m×m symmetric tridiagonal matrix represented by the diagonal d
    75  // and subdiagonal e, and B and C are m×n matrices.
    76  func dstmm(m, n int, d, e []float64, b []float64, ldb int, c []float64, ldc int) {
    77  	if m == 0 || n == 0 {
    78  		return
    79  	}
    80  	if m == 1 {
    81  		d0 := d[0]
    82  		for j, b0j := range b[:n] {
    83  			c[j] = d0 * b0j
    84  		}
    85  		return
    86  	}
    87  	for j := 0; j < n; j++ {
    88  		c[j] = d[0]*b[j] + e[0]*b[ldb+j]
    89  	}
    90  	for i := 1; i < m-1; i++ {
    91  		for j := 0; j < n; j++ {
    92  			c[i*ldc+j] = e[i-1]*b[(i-1)*ldb+j] + d[i]*b[i*ldb+j] + e[i]*b[(i+1)*ldb+j]
    93  		}
    94  	}
    95  	for j := 0; j < n; j++ {
    96  		c[(m-1)*ldc+j] = e[m-2]*b[(m-2)*ldb+j] + d[m-1]*b[(m-1)*ldb+j]
    97  	}
    98  }
    99  
   100  // dpttrsResidual returns |XGOT - XWANT|_1 / n.
   101  func dpttrsResidual(xGot, xWant blas64.General) float64 {
   102  	n, nrhs := xGot.Rows, xGot.Cols
   103  	d := zeros(n, nrhs, nrhs)
   104  	for i := 0; i < n; i++ {
   105  		for j := 0; j < nrhs; j++ {
   106  			d.Data[i*d.Stride+j] = xGot.Data[i*xGot.Stride+j] - xWant.Data[i*xWant.Stride+j]
   107  		}
   108  	}
   109  	return dlange(lapack.MaxColumnSum, n, nrhs, d.Data, d.Stride) / float64(n)
   110  }