+++ /dev/null
-// 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
-}
OSDN Git Service
