gonum.org/v1/gonum@v0.14.0/lapack/testlapack/dlatrs.go (about)

     1  // Copyright ©2017 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/blas"
    15  	"gonum.org/v1/gonum/blas/blas64"
    16  )
    17  
    18  type Dlatrser interface {
    19  	Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, normin bool, n int, a []float64, lda int, x []float64, cnorm []float64) (scale float64)
    20  }
    21  
    22  func DlatrsTest(t *testing.T, impl Dlatrser) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    25  		for _, trans := range []blas.Transpose{blas.Trans, blas.NoTrans} {
    26  			for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 7, 10, 20, 50, 100} {
    27  				for _, lda := range []int{n, 2*n + 1} {
    28  					lda = max(1, lda)
    29  					imats := []int{7, 11, 12, 13, 14, 15, 16, 17, 18, 19}
    30  					if n < 6 {
    31  						imats = append(imats, 19)
    32  					}
    33  					for _, imat := range imats {
    34  						testDlatrs(t, impl, imat, uplo, trans, n, lda, rnd)
    35  					}
    36  				}
    37  			}
    38  		}
    39  	}
    40  }
    41  
    42  func testDlatrs(t *testing.T, impl Dlatrser, imat int, uplo blas.Uplo, trans blas.Transpose, n, lda int, rnd *rand.Rand) {
    43  	const tol = 1e-14
    44  
    45  	a := nanSlice(n * lda)
    46  	b := nanSlice(n)
    47  	work := make([]float64, 3*n)
    48  
    49  	// Generate triangular test matrix and right hand side.
    50  	diag := dlattr(imat, uplo, trans, n, a, lda, b, work, rnd)
    51  	if imat <= 10 {
    52  		// b has not been generated.
    53  		dlarnv(b, 3, rnd)
    54  	}
    55  
    56  	cnorm := nanSlice(n)
    57  	x := make([]float64, n)
    58  
    59  	// Call Dlatrs with normin=false.
    60  	copy(x, b)
    61  	scale := impl.Dlatrs(uplo, trans, diag, false, n, a, lda, x, cnorm)
    62  	prefix := fmt.Sprintf("Case imat=%v (n=%v,lda=%v,trans=%c,uplo=%c,diag=%c", imat, n, lda, trans, uplo, diag)
    63  	for i, v := range cnorm {
    64  		if math.IsNaN(v) {
    65  			t.Errorf("%v: cnorm[%v] not computed (scale=%v,normin=false)", prefix, i, scale)
    66  		}
    67  	}
    68  	resid, hasNaN := dlatrsResidual(uplo, trans, diag, n, a, lda, scale, cnorm, x, b, work[:n])
    69  	if hasNaN {
    70  		t.Errorf("%v: unexpected NaN (scale=%v,normin=false)", prefix, scale)
    71  	} else if resid > tol {
    72  		t.Errorf("%v: residual %v too large (scale=%v,normin=false)", prefix, resid, scale)
    73  	}
    74  
    75  	// Call Dlatrs with normin=true because cnorm has been filled.
    76  	copy(x, b)
    77  	scale = impl.Dlatrs(uplo, trans, diag, true, n, a, lda, x, cnorm)
    78  	resid, hasNaN = dlatrsResidual(uplo, trans, diag, n, a, lda, scale, cnorm, x, b, work[:n])
    79  	if hasNaN {
    80  		t.Errorf("%v: unexpected NaN (scale=%v,normin=true)", prefix, scale)
    81  	} else if resid > tol {
    82  		t.Errorf("%v: residual %v too large (scale=%v,normin=true)", prefix, resid, scale)
    83  	}
    84  }
    85  
    86  // dlatrsResidual returns norm(trans(A)*x-scale*b) / (norm(trans(A))*norm(x)*eps)
    87  // and whether NaN has been encountered in the process.
    88  func dlatrsResidual(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []float64, lda int, scale float64, cnorm []float64, x, b, work []float64) (resid float64, hasNaN bool) {
    89  	if n == 0 {
    90  		return 0, false
    91  	}
    92  
    93  	// Compute the norm of the triangular matrix A using the column norms
    94  	// already computed by Dlatrs.
    95  	var tnorm float64
    96  	if diag == blas.NonUnit {
    97  		for j := 0; j < n; j++ {
    98  			tnorm = math.Max(tnorm, math.Abs(a[j*lda+j])+cnorm[j])
    99  		}
   100  	} else {
   101  		for j := 0; j < n; j++ {
   102  			tnorm = math.Max(tnorm, 1+cnorm[j])
   103  		}
   104  	}
   105  
   106  	const (
   107  		eps  = dlamchE
   108  		tiny = safmin
   109  	)
   110  
   111  	bi := blas64.Implementation()
   112  
   113  	// Compute norm(trans(A)*x-scale*b) / (norm(trans(A))*norm(x)*eps)
   114  	copy(work, x)
   115  	ix := bi.Idamax(n, work, 1)
   116  	xnorm := math.Max(1, math.Abs(work[ix]))
   117  	xscal := 1 / xnorm / float64(n)
   118  	bi.Dscal(n, xscal, work, 1)
   119  	bi.Dtrmv(uplo, trans, diag, n, a, lda, work, 1)
   120  	bi.Daxpy(n, -scale*xscal, b, 1, work, 1)
   121  	for _, v := range work {
   122  		if math.IsNaN(v) {
   123  			return 1 / eps, true
   124  		}
   125  	}
   126  	ix = bi.Idamax(n, work, 1)
   127  	resid = math.Abs(work[ix])
   128  	ix = bi.Idamax(n, x, 1)
   129  	xnorm = math.Abs(x[ix])
   130  	if resid*tiny <= xnorm {
   131  		if xnorm > 0 {
   132  			resid /= xnorm
   133  		}
   134  	} else if resid > 0 {
   135  		resid = 1 / eps
   136  	}
   137  	if resid*tiny <= tnorm {
   138  		if tnorm > 0 {
   139  			resid /= tnorm
   140  		}
   141  	} else if resid > 0 {
   142  		resid = 1 / eps
   143  	}
   144  	return resid, false
   145  }