testlapack: consolidate floating-point constants
MODIFY
lapack/testlapack/dgeev.go
lapack/testlapack/dgeev.go
@@ -769,10 +769,9 @@ func residualRightEV(a, e blas64.General, wr, wi []float64) float64 {
}
bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, a.Data, a.Stride, e.Data, e.Stride, -1, r, ldr)
- unfl := dlamchS
- ulp := dlamchE
- anorm := math.Max(dlange(lapack.MaxColumnSum, n, n, a.Data, a.Stride), unfl)
- enorm := math.Max(dlange(lapack.MaxColumnSum, n, n, e.Data, e.Stride), ulp)
+ const eps = dlamchE
+ anorm := math.Max(dlange(lapack.MaxColumnSum, n, n, a.Data, a.Stride), safmin)
+ enorm := math.Max(dlange(lapack.MaxColumnSum, n, n, e.Data, e.Stride), eps)
errnorm := dlange(lapack.MaxColumnSum, n, n, r, ldr) / enorm
if anorm > errnorm {
return errnorm / anorm
@@ -824,10 +823,9 @@ func residualLeftEV(a, e blas64.General, wr, wi []float64) float64 {
}
bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, a.Data, a.Stride, e.Data, e.Stride, -1, r, ldr)
- unfl := dlamchS
- ulp := dlamchE
- anorm := math.Max(dlange(lapack.MaxRowSum, n, n, a.Data, a.Stride), unfl)
- enorm := math.Max(dlange(lapack.MaxColumnSum, n, n, e.Data, e.Stride), ulp)
+ const eps = dlamchE
+ anorm := math.Max(dlange(lapack.MaxRowSum, n, n, a.Data, a.Stride), safmin)
+ enorm := math.Max(dlange(lapack.MaxColumnSum, n, n, e.Data, e.Stride), eps)
errnorm := dlange(lapack.MaxColumnSum, n, n, r, ldr) / enorm
if anorm > errnorm {
return errnorm / anorm

MODIFY
lapack/testlapack/dgesvd.go
lapack/testlapack/dgesvd.go
@@ -96,15 +96,12 @@ func dgesvdTest(t *testing.T, impl Dgesvder, m, n, mtype int, tol float64) {
false, // signs of s[i] are not randomly flipped
1, rnd) // random numbers are drawn uniformly from [0,1)
// Decide scale factor for the singular values based on the matrix type.
- ulp := dlamchP
- unfl := dlamchS
- ovfl := 1 / unfl
aNorm = 1
if mtype == 4 {
- aNorm = unfl / ulp
+ aNorm = smlnum
}
if mtype == 5 {
- aNorm = ovfl * ulp
+ aNorm = bignum
}
// Scale singular values so that the maximum singular value is
// equal to aNorm (we know that the singular values are

MODIFY
lapack/testlapack/dlangb_bench.go
lapack/testlapack/dlangb_bench.go
@@ -14,14 +14,6 @@ import (
)
func DlangbBenchmark(b *testing.B, impl Dlangber) {
- const (
- safmin = dlamchS
- safmax = 1 / safmin
- ulp = dlamchP
- smlnum = safmin / ulp
- bignum = safmax * ulp
- )
-
rnd := rand.New(rand.NewSource(1))
for _, bm := range []struct {
n, k int

MODIFY
lapack/testlapack/dlartg.go
lapack/testlapack/dlartg.go
@@ -15,23 +15,21 @@ type Dlartger interface {
}
func DlartgTest(t *testing.T, impl Dlartger) {
- const tol = 20 * dlamchP
+ const tol = 20 * ulp
- safmin := dlamchS
- safmax := 1 / safmin
values := []float64{
-safmax,
- -1 / dlamchP,
+ -1 / ulp,
-1,
-1.0 / 3,
- -dlamchP,
+ -ulp,
-safmin,
0,
safmin,
- dlamchP,
+ ulp,
1.0 / 3,
1,
- 1 / dlamchP,
+ 1 / ulp,
safmax,
math.Inf(-1),
math.Inf(1),

MODIFY
lapack/testlapack/dlassq.go
lapack/testlapack/dlassq.go
@@ -19,13 +19,6 @@ type Dlassqer interface {
}
func DlassqTest(t *testing.T, impl Dlassqer) {
- const (
- safmin = dlamchS
- safmax = 1 / safmin
- ulp = dlamchP
- smlnum = safmin / ulp
- bignum = safmax * ulp
- )
values := []float64{
0,
2 * safmin,

MODIFY
lapack/testlapack/dlatbs.go
lapack/testlapack/dlatbs.go
@@ -134,9 +134,12 @@ func dlatbsResidual(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, kd
}
}
+ const (
+ eps = dlamchE
+ tiny = safmin
+ )
+
bi := blas64.Implementation()
- eps := dlamchE
- smlnum := dlamchS
ix := bi.Idamax(n, x, 1)
xNorm := math.Max(1, math.Abs(x[ix]))
@@ -150,14 +153,14 @@ func dlatbsResidual(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, kd
ix = bi.Idamax(n, resid, 1)
residNorm := math.Abs(resid[ix])
- if residNorm*smlnum <= xNorm {
+ if residNorm*tiny <= xNorm {
if xNorm > 0 {
residNorm /= xNorm
}
} else if residNorm > 0 {
residNorm = 1 / eps
}
- if residNorm*smlnum <= tnorm {
+ if residNorm*tiny <= tnorm {
if tnorm > 0 {
residNorm /= tnorm
}

MODIFY
lapack/testlapack/dlatrs.go
lapack/testlapack/dlatrs.go
@@ -103,8 +103,11 @@ func dlatrsResidual(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int,
}
}
- eps := dlamchE
- smlnum := dlamchS
+ const (
+ eps = dlamchE
+ tiny = safmin
+ )
+
bi := blas64.Implementation()
// Compute norm(trans(A)*x-scale*b) / (norm(trans(A))*norm(x)*eps)
@@ -124,14 +127,14 @@ func dlatrsResidual(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int,
resid = math.Abs(work[ix])
ix = bi.Idamax(n, x, 1)
xnorm = math.Abs(x[ix])
- if resid*smlnum <= xnorm {
+ if resid*tiny <= xnorm {
if xnorm > 0 {
resid /= xnorm
}
} else if resid > 0 {
resid = 1 / eps
}
- if resid*smlnum <= tnorm {
+ if resid*tiny <= tnorm {
if tnorm > 0 {
resid /= tnorm
}

MODIFY
lapack/testlapack/general.go
lapack/testlapack/general.go
@@ -24,6 +24,12 @@ const (
dlamchP = dlamchB * dlamchE
// dlamchS is the smallest normal number. For IEEE this is 2^{-1022}.
dlamchS = 0x1p-1022
+
+ safmin = dlamchS
+ safmax = 1 / safmin
+ ulp = dlamchP
+ smlnum = safmin / ulp
+ bignum = safmax * ulp
)
func max(a, b int) int {

MODIFY
lapack/testlapack/matgen.go
lapack/testlapack/matgen.go
@@ -341,9 +341,10 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
return blas.NonUnit
}
- ulp := dlamchE * dlamchB
- smlnum := dlamchS
- bignum := (1 - ulp) / smlnum
+ const (
+ tiny = safmin
+ huge = (1 - ulp) / tiny
+ )
bi := blas64.Implementation()
@@ -388,10 +389,10 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
for i := 0; i < n; i++ {
a[i*lda+i] = math.Copysign(2, a[i*lda+i])
}
- // Set the right hand side so that the largest value is bignum.
+ // Set the right hand side so that the largest value is huge.
dlarnv(b, 2, rnd)
imax := bi.Idamax(n, b, 1)
- bscal := bignum / math.Max(1, b[imax])
+ bscal := huge / math.Max(1, b[imax])
bi.Dscal(n, bscal, b, 1)
case 12:
// Make the first diagonal element in the solve small to cause
@@ -406,14 +407,14 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
bi.Dscal(n-i-1, tscal, a[i*lda+i+1:], 1)
a[i*lda+i] = math.Copysign(1, a[i*lda+i])
}
- a[(n-1)*lda+n-1] *= smlnum
+ a[(n-1)*lda+n-1] *= tiny
case blas.Lower:
for i := 0; i < n; i++ {
dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
bi.Dscal(i, tscal, a[i*lda:], 1)
a[i*lda+i] = math.Copysign(1, a[i*lda+i])
}
- a[0] *= smlnum
+ a[0] *= tiny
}
dlarnv(b, 2, rnd)
case 13:
@@ -427,13 +428,13 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
dlarnv(a[i*lda+i:i*lda+n], 2, rnd)
a[i*lda+i] = math.Copysign(1, a[i*lda+i])
}
- a[(n-1)*lda+n-1] *= smlnum
+ a[(n-1)*lda+n-1] *= tiny
case blas.Lower:
for i := 0; i < n; i++ {
dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
a[i*lda+i] = math.Copysign(1, a[i*lda+i])
}
- a[0] *= smlnum
+ a[0] *= tiny
}
dlarnv(b, 2, rnd)
case 14:
@@ -448,7 +449,7 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
a[i*lda+j] = 0
}
if (n-1-i)&0x2 == 0 {
- a[i*lda+i] = smlnum
+ a[i*lda+i] = tiny
} else {
a[i*lda+i] = 1
}
@@ -459,7 +460,7 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
a[i*lda+j] = 0
}
if i&0x2 == 0 {
- a[i*lda+i] = smlnum
+ a[i*lda+i] = tiny
} else {
a[i*lda+i] = 1
}
@@ -471,12 +472,12 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
b[0] = 0
for i := n - 1; i > 0; i -= 2 {
b[i] = 0
- b[i-1] = smlnum
+ b[i-1] = tiny
}
case blas.Lower:
for i := 0; i < n-1; i += 2 {
b[i] = 0
- b[i+1] = smlnum
+ b[i+1] = tiny
}
b[n-1] = 0
}
@@ -486,7 +487,7 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
// needed, the matrix is bidiagonal.
diag = blas.NonUnit
texp := 1 / math.Max(1, float64(n-1))
- tscal := math.Pow(smlnum, texp)
+ tscal := math.Pow(tiny, texp)
switch uplo {
case blas.Upper:
for i := 0; i < n; i++ {
@@ -535,7 +536,7 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
// is constructed to cause overflow when adding a column in
// every other step.
diag = blas.NonUnit
- tscal := (1 - ulp) / (dlamchS / ulp)
+ tscal := (1 - ulp) / tiny
texp := 1.0
switch uplo {
case blas.Upper:
@@ -588,20 +589,20 @@ func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64,
a[i*lda+i] = math.NaN()
}
}
- // Set the right hand side so that the largest value is bignum.
+ // Set the right hand side so that the largest value is huge.
dlarnv(b, 2, rnd)
iy := bi.Idamax(n, b, 1)
bnorm := math.Abs(b[iy])
- bscal := bignum / math.Max(1, bnorm)
+ bscal := huge / math.Max(1, bnorm)
bi.Dscal(n, bscal, b, 1)
case 19:
// Generate a triangular matrix with elements between
- // bignum/(n-1) and bignum so that at least one of the column
- // norms will exceed bignum.
+ // huge/(n-1) and huge so that at least one of the column
+ // norms will exceed huge.
// Dlatrs cannot handle this case for (typically) n>5.
diag = blas.NonUnit
- tleft := bignum / math.Max(1, float64(n-1))
- tscal := bignum * (float64(n-1) / math.Max(1, float64(n)))
+ tleft := huge / math.Max(1, float64(n-1))
+ tscal := huge * (float64(n-1) / math.Max(1, float64(n)))
switch uplo {
case blas.Upper:
for i := 0; i < n; i++ {
@@ -754,15 +755,12 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
}
const (
- unfl = dlamchS
- ulp = dlamchE * dlamchB
- smlnum = unfl
- bignum = (1 - ulp) / smlnum
+ tiny = safmin
+ huge = (1 - ulp) / tiny
- eps = dlamchP
- small = 0.25 * (dlamchS / eps)
+ small = 0.25 * (safmin / ulp)
large = 1 / small
- badc2 = 0.1 / eps
+ badc2 = 0.1 / ulp
)
badc1 := math.Sqrt(badc2)
@@ -969,9 +967,9 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
ab[i*ldab+kd] = math.Copysign(2, ab[i*ldab+kd])
}
}
- // Set the right hand side so that the largest value is bignum.
+ // Set the right hand side so that the largest value is huge.
bnorm := math.Abs(b[bi.Idamax(n, b, 1)])
- bscal := bignum / math.Max(1, bnorm)
+ bscal := huge / math.Max(1, bnorm)
bi.Dscal(n, bscal, b, 1)
case 11:
@@ -989,7 +987,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
}
ab[i*ldab] = math.Copysign(1, ab[i*ldab])
}
- ab[(n-1)*ldab] *= smlnum
+ ab[(n-1)*ldab] *= tiny
} else {
for i := 0; i < n; i++ {
jlen := min(i+1, kd+1)
@@ -1000,7 +998,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
}
ab[i*ldab+kd] = math.Copysign(1, ab[i*ldab+kd])
}
- ab[kd] *= smlnum
+ ab[kd] *= tiny
}
case 12:
@@ -1014,7 +1012,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
dlarnv(arow, 2, rnd)
ab[i*ldab] = math.Copysign(1, ab[i*ldab])
}
- ab[(n-1)*ldab] *= smlnum
+ ab[(n-1)*ldab] *= tiny
} else {
for i := 0; i < n; i++ {
jlen := min(i+1, kd+1)
@@ -1022,7 +1020,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
dlarnv(arow, 2, rnd)
ab[i*ldab+kd] = math.Copysign(1, ab[i*ldab+kd])
}
- ab[kd] *= smlnum
+ ab[kd] *= tiny
}
case 13:
@@ -1033,7 +1031,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
icount := 1
for i := n - 1; i >= 0; i-- {
if icount <= 2 {
- ab[i*ldab] = smlnum
+ ab[i*ldab] = tiny
} else {
ab[i*ldab] = 1
}
@@ -1052,7 +1050,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
ab[i*ldab+j] = 0
}
if icount <= 2 {
- ab[i*ldab+kd] = smlnum
+ ab[i*ldab+kd] = tiny
} else {
ab[i*ldab+kd] = 1
}
@@ -1067,13 +1065,13 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
b[0] = 0
for i := n - 1; i > 1; i -= 2 {
b[i] = 0
- b[i-1] = smlnum
+ b[i-1] = tiny
}
} else {
b[n-1] = 0
for i := 0; i < n-1; i += 2 {
b[i] = 0
- b[i+1] = smlnum
+ b[i+1] = tiny
}
}
@@ -1081,7 +1079,7 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
// Make the diagonal elements small to cause gradual overflow when
// dividing by T[j,j]. To control the amount of scaling needed, the
// matrix is bidiagonal.
- tscal := math.Pow(smlnum, 1/float64(kd+1))
+ tscal := math.Pow(tiny, 1/float64(kd+1))
if uplo == blas.Upper {
for i := 0; i < n; i++ {
ab[i*ldab] = tscal
@@ -1191,21 +1189,21 @@ func dlattb(kind int, uplo blas.Uplo, trans blas.Transpose, n, kd int, ab []floa
ab[i*ldab+kd] = float64(i + 2)
}
}
- // Set the right hand side so that the largest value is bignum.
+ // Set the right hand side so that the largest value is huge.
bnorm := math.Abs(b[bi.Idamax(n, b, 1)])
- bscal := bignum / math.Max(1, bnorm)
+ bscal := huge / math.Max(1, bnorm)
bi.Dscal(n, bscal, b, 1)
case 18:
- // Generate a triangular matrix with elements between bignum/kd and
- // bignum so that at least one of the column norms will exceed bignum.
- tleft := bignum / math.Max(1, float64(kd))
+ // Generate a triangular matrix with elements between huge/kd and
+ // huge so that at least one of the column norms will exceed huge.
+ tleft := huge / math.Max(1, float64(kd))
// The reference LAPACK has
- // tscal := bignum * (float64(kd) / float64(kd+1))
+ // tscal := huge * (float64(kd) / float64(kd+1))
// but this causes overflow when computing cnorm in Dlatbs. Our choice
// is more conservative but increases coverage in the same way as the
// LAPACK version.
- tscal := bignum / math.Max(1, float64(kd))
+ tscal := huge / math.Max(1, float64(kd))
if uplo == blas.Upper {
for i := 0; i < n; i++ {
for j := 0; j < min(n-i, kd+1); j++ {

testlapack: consolidate floating-point constants