1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
// 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, 19}
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=%c,uplo=%c,diag=%c", 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])
}
}
const (
eps = dlamchE
tiny = safmin
)
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*tiny <= xnorm {
if xnorm > 0 {
resid /= xnorm
}
} else if resid > 0 {
resid = 1 / eps
}
if resid*tiny <= tnorm {
if tnorm > 0 {
resid /= tnorm
}
} else if resid > 0 {
resid = 1 / eps
}
return resid, false
}
|
