Hulk did something

[bytom/vapor.git] / vendor / gonum.org / v1 / gonum / lapack / testlapack / dlatrs.go

diff --git a/vendor/gonum.org/v1/gonum/lapack/testlapack/dlatrs.go b/vendor/gonum.org/v1/gonum/lapack/testlapack/dlatrs.go
new file mode 100644 (file)
index 0000000..6d0bbda
--- /dev/null
+++ b/vendor/gonum.org/v1/gonum/lapack/testlapack/dlatrs.go
@@ -0,0 +1,142 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package testlapack
+
+import (
+       "fmt"
+       "math"
+       "testing"
+
+       "golang.org/x/exp/rand"
+
+       "gonum.org/v1/gonum/blas"
+       "gonum.org/v1/gonum/blas/blas64"
+)
+
+type Dlatrser interface {
+       Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, normin bool, n int, a []float64, lda int, x []float64, cnorm []float64) (scale float64)
+}
+
+func DlatrsTest(t *testing.T, impl Dlatrser) {
+       rnd := rand.New(rand.NewSource(1))
+       for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
+               for _, trans := range []blas.Transpose{blas.Trans, blas.NoTrans} {
+                       for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 7, 10, 20, 50, 100} {
+                               for _, lda := range []int{n, 2*n + 1} {
+                                       lda = max(1, lda)
+                                       imats := []int{7, 11, 12, 13, 14, 15, 16, 17, 18}
+                                       if n < 6 {
+                                               imats = append(imats, 19)
+                                       }
+                                       for _, imat := range imats {
+                                               testDlatrs(t, impl, imat, uplo, trans, n, lda, rnd)
+                                       }
+                               }
+                       }
+               }
+       }
+}
+
+func testDlatrs(t *testing.T, impl Dlatrser, imat int, uplo blas.Uplo, trans blas.Transpose, n, lda int, rnd *rand.Rand) {
+       const tol = 1e-14
+
+       a := nanSlice(n * lda)
+       b := nanSlice(n)
+       work := make([]float64, 3*n)
+
+       // Generate triangular test matrix and right hand side.
+       diag := dlattr(imat, uplo, trans, n, a, lda, b, work, rnd)
+       if imat <= 10 {
+               // b has not been generated.
+               dlarnv(b, 3, rnd)
+       }
+
+       cnorm := nanSlice(n)
+       x := make([]float64, n)
+
+       // Call Dlatrs with normin=false.
+       copy(x, b)
+       scale := impl.Dlatrs(uplo, trans, diag, false, n, a, lda, x, cnorm)
+       prefix := fmt.Sprintf("Case imat=%v (n=%v,lda=%v,trans=%v,uplo=%v,diag=%v", imat, n, lda, trans, uplo, diag)
+       for i, v := range cnorm {
+               if math.IsNaN(v) {
+                       t.Errorf("%v: cnorm[%v] not computed (scale=%v,normin=false)", prefix, i, scale)
+               }
+       }
+       resid, hasNaN := dlatrsResidual(uplo, trans, diag, n, a, lda, scale, cnorm, x, b, work[:n])
+       if hasNaN {
+               t.Errorf("%v: unexpected NaN (scale=%v,normin=false)", prefix, scale)
+       } else if resid > tol {
+               t.Errorf("%v: residual %v too large (scale=%v,normin=false)", prefix, resid, scale)
+       }
+
+       // Call Dlatrs with normin=true because cnorm has been filled.
+       copy(x, b)
+       scale = impl.Dlatrs(uplo, trans, diag, true, n, a, lda, x, cnorm)
+       resid, hasNaN = dlatrsResidual(uplo, trans, diag, n, a, lda, scale, cnorm, x, b, work[:n])
+       if hasNaN {
+               t.Errorf("%v: unexpected NaN (scale=%v,normin=true)", prefix, scale)
+       } else if resid > tol {
+               t.Errorf("%v: residual %v too large (scale=%v,normin=true)", prefix, resid, scale)
+       }
+}
+
+// dlatrsResidual returns norm(trans(A)*x-scale*b) / (norm(trans(A))*norm(x)*eps)
+// and whether NaN has been encountered in the process.
+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) {
+       if n == 0 {
+               return 0, false
+       }
+
+       // Compute the norm of the triangular matrix A using the column norms
+       // already computed by Dlatrs.
+       var tnorm float64
+       if diag == blas.NonUnit {
+               for j := 0; j < n; j++ {
+                       tnorm = math.Max(tnorm, math.Abs(a[j*lda+j])+cnorm[j])
+               }
+       } else {
+               for j := 0; j < n; j++ {
+                       tnorm = math.Max(tnorm, 1+cnorm[j])
+               }
+       }
+
+       eps := dlamchE
+       smlnum := dlamchS
+       bi := blas64.Implementation()
+
+       // Compute norm(trans(A)*x-scale*b) / (norm(trans(A))*norm(x)*eps)
+       copy(work, x)
+       ix := bi.Idamax(n, work, 1)
+       xnorm := math.Max(1, math.Abs(work[ix]))
+       xscal := 1 / xnorm / float64(n)
+       bi.Dscal(n, xscal, work, 1)
+       bi.Dtrmv(uplo, trans, diag, n, a, lda, work, 1)
+       bi.Daxpy(n, -scale*xscal, b, 1, work, 1)
+       for _, v := range work {
+               if math.IsNaN(v) {
+                       return 1 / eps, true
+               }
+       }
+       ix = bi.Idamax(n, work, 1)
+       resid = math.Abs(work[ix])
+       ix = bi.Idamax(n, x, 1)
+       xnorm = math.Abs(x[ix])
+       if resid*smlnum <= xnorm {
+               if xnorm > 0 {
+                       resid /= xnorm
+               }
+       } else if resid > 0 {
+               resid = 1 / eps
+       }
+       if resid*smlnum <= tnorm {
+               if tnorm > 0 {
+                       resid /= tnorm
+               }
+       } else if resid > 0 {
+               resid = 1 / eps
+       }
+       return resid, false
+}