github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/testlapack/dtrtrs.go (about)

     1  // Copyright ©2020 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  	"github.com/jingcheng-WU/gonum/blas"
    15  	"github.com/jingcheng-WU/gonum/blas/blas64"
    16  	"github.com/jingcheng-WU/gonum/floats"
    17  	"github.com/jingcheng-WU/gonum/lapack"
    18  )
    19  
    20  type Dtrtrser interface {
    21  	Dtrtrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, a []float64, lda int, b []float64, ldb int) bool
    22  }
    23  
    24  func DtrtrsTest(t *testing.T, impl Dtrtrser) {
    25  	rnd := rand.New(rand.NewSource(1))
    26  
    27  	for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans, blas.ConjTrans} {
    28  		name := transToString(trans)
    29  		t.Run(name, func(t *testing.T) {
    30  			for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    31  				for _, diag := range []blas.Diag{blas.Unit, blas.NonUnit} {
    32  					for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
    33  						for _, nrhs := range []int{0, 1, 2, 3, 4, 5, 10, 15} {
    34  							for _, lda := range []int{max(1, n), n + 3} {
    35  								for _, ldb := range []int{max(1, nrhs), nrhs + 3} {
    36  									if diag == blas.Unit {
    37  										dtrtrsTest(t, impl, rnd, uplo, trans, diag, n, nrhs, lda, ldb, false)
    38  									} else {
    39  										dtrtrsTest(t, impl, rnd, uplo, trans, diag, n, nrhs, lda, ldb, true)
    40  										dtrtrsTest(t, impl, rnd, uplo, trans, diag, n, nrhs, lda, ldb, false)
    41  									}
    42  								}
    43  							}
    44  						}
    45  					}
    46  				}
    47  			}
    48  		})
    49  	}
    50  }
    51  
    52  func dtrtrsTest(t *testing.T, impl Dtrtrser, rnd *rand.Rand, uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, lda, ldb int, singular bool) {
    53  	if singular && diag == blas.Unit {
    54  		panic("blas.Unit triangular matrix cannot be singular")
    55  	}
    56  
    57  	const tol = 1e-14
    58  
    59  	if n == 0 {
    60  		singular = false
    61  	}
    62  	name := fmt.Sprintf("uplo=%v,diag=%v,n=%v,nrhs=%v,lda=%v,ldb=%v,sing=%v", string(uplo), string(diag), n, nrhs, lda, ldb, singular)
    63  
    64  	// Generate a random triangular matrix A. One of its triangles won't be
    65  	// referenced.
    66  	a := make([]float64, n*lda)
    67  	for i := range a {
    68  		a[i] = rnd.NormFloat64()
    69  	}
    70  	if singular {
    71  		i := rnd.Intn(n)
    72  		a[i*lda+i] = 0
    73  	}
    74  	aCopy := make([]float64, len(a))
    75  	copy(aCopy, a)
    76  
    77  	// Generate a random solution matrix X.
    78  	x := make([]float64, n*ldb)
    79  	for i := range x {
    80  		x[i] = rnd.NormFloat64()
    81  	}
    82  
    83  	// Generate the right-hand side as A * X  or  Aᵀ * X.
    84  	b := make([]float64, len(x))
    85  	copy(b, x)
    86  	bi := blas64.Implementation()
    87  	bi.Dtrmm(blas.Left, uplo, trans, diag, n, nrhs, 1, a, lda, b, ldb)
    88  
    89  	got := make([]float64, len(b))
    90  	copy(got, b)
    91  
    92  	ok := impl.Dtrtrs(uplo, trans, diag, n, nrhs, a, lda, got, ldb)
    93  
    94  	if !floats.Equal(a, aCopy) {
    95  		t.Errorf("%v: unexpected modification of A", name)
    96  	}
    97  
    98  	if ok == singular {
    99  		t.Errorf("%v: misdetected singular matrix, ok=%v", name, ok)
   100  	}
   101  
   102  	if !ok {
   103  		if !floats.Equal(got, b) {
   104  			t.Errorf("%v: unexpected modification of B when singular", name)
   105  		}
   106  		return
   107  	}
   108  
   109  	if n == 0 || nrhs == 0 {
   110  		return
   111  	}
   112  
   113  	work := make([]float64, n)
   114  
   115  	// Compute the 1-norm of A or Aᵀ.
   116  	var aNorm float64
   117  	if trans == blas.NoTrans {
   118  		aNorm = dlantr(lapack.MaxColumnSum, uplo, diag, n, n, a, lda, work)
   119  	} else {
   120  		aNorm = dlantr(lapack.MaxRowSum, uplo, diag, n, n, a, lda, work)
   121  	}
   122  
   123  	// Compute the maximum over the number of right-hand sides of
   124  	//  |op(A)*x-b| / (|op(A)| * |x|)
   125  	var resid float64
   126  	for j := 0; j < nrhs; j++ {
   127  		bi.Dcopy(n, got[j:], ldb, work, 1)
   128  		bi.Dtrmv(uplo, trans, diag, n, a, lda, work, 1)
   129  		bi.Daxpy(n, -1, b[j:], ldb, work, 1)
   130  		rjNorm := bi.Dasum(n, work, 1)
   131  		xNorm := bi.Dasum(n, got[j:], ldb)
   132  		resid = math.Max(resid, rjNorm/aNorm/xNorm)
   133  	}
   134  	if resid > tol {
   135  		t.Errorf("%v: unexpected result; resid=%v,want<=%v", name, resid, tol)
   136  	}
   137  }