vendor/gonum.org/v1/gonum/lapack/testlapack/dlatrs.go

   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}
  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=%v,uplo=%v,diag=%v", 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         eps := dlamchE
 107         smlnum := dlamchS
 108         bi := blas64.Implementation()
 109
 110         // Compute norm(trans(A)*x-scale*b) / (norm(trans(A))*norm(x)*eps)
 111         copy(work, x)
 112         ix := bi.Idamax(n, work, 1)
 113         xnorm := math.Max(1, math.Abs(work[ix]))
 114         xscal := 1 / xnorm / float64(n)
 115         bi.Dscal(n, xscal, work, 1)
 116         bi.Dtrmv(uplo, trans, diag, n, a, lda, work, 1)
 117         bi.Daxpy(n, -scale*xscal, b, 1, work, 1)
 118         for _, v := range work {
 119                 if math.IsNaN(v) {
 120                         return 1 / eps, true
 121                 }
 122         }
 123         ix = bi.Idamax(n, work, 1)
 124         resid = math.Abs(work[ix])
 125         ix = bi.Idamax(n, x, 1)
 126         xnorm = math.Abs(x[ix])
 127         if resid*smlnum <= xnorm {
 128                 if xnorm > 0 {
 129                         resid /= xnorm
 130                 }
 131         } else if resid > 0 {
 132                 resid = 1 / eps
 133         }
 134         if resid*smlnum <= tnorm {
 135                 if tnorm > 0 {
 136                         resid /= tnorm
 137                 }
 138         } else if resid > 0 {
 139                 resid = 1 / eps
 140         }
 141         return resid, false
 142 }