lapack/testlapack: extend test for Dlasy2 to cover rescaling
MODIFY
lapack/testlapack/dlasy2.go
lapack/testlapack/dlasy2.go
@@ -6,13 +6,13 @@ package testlapack
import (
"fmt"
- "math"
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
+ "gonum.org/v1/gonum/lapack"
)
type Dlasy2er interface {
@@ -26,10 +26,13 @@ func Dlasy2Test(t *testing.T, impl Dlasy2er) {
for _, isgn := range []int{1, -1} {
for _, n1 := range []int{0, 1, 2} {
for _, n2 := range []int{0, 1, 2} {
- for _, extra := range []int{0, 1, 2, 13} {
- for cas := 0; cas < 1000; cas++ {
- testDlasy2(t, impl, tranl, tranr, isgn, n1, n2, extra, rnd)
+ for _, extra := range []int{0, 3} {
+ for cas := 0; cas < 100; cas++ {
+ var big bool
+ if cas%2 == 0 {
+ big = true
}
+ testDlasy2(t, impl, tranl, tranr, isgn, n1, n2, extra, big, rnd)
}
}
}
@@ -37,69 +40,85 @@ func Dlasy2Test(t *testing.T, impl Dlasy2er) {
}
}
}
+}
+
+func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, big bool, rnd *rand.Rand) {
+ const tol = 1e-14
-func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, rnd *rand.Rand) {
- const tol = 1e-10
+ name := fmt.Sprintf("Case n1=%v, n2=%v, isgn=%v, big=%v", n1, n2, isgn, big)
tl := randomGeneral(n1, n1, n1+extra, rnd)
tr := randomGeneral(n2, n2, n2+extra, rnd)
- b := randomGeneral(n1, n2, n2+extra, rnd)
x := randomGeneral(n1, n2, n2+extra, rnd)
+ b := randomGeneral(n1, n2, n2+extra, rnd)
+ if big {
+ for i := 0; i < n1; i++ {
+ for j := 0; j < n2; j++ {
+ b.Data[i*b.Stride+j] *= bignum
+ }
+ }
+ }
+
+ tlCopy := cloneGeneral(tl)
+ trCopy := cloneGeneral(tr)
+ bCopy := cloneGeneral(b)
scale, xnorm, ok := impl.Dlasy2(tranl, tranr, isgn, n1, n2, tl.Data, tl.Stride, tr.Data, tr.Stride, b.Data, b.Stride, x.Data, x.Stride)
- if scale > 1 {
- t.Errorf("invalid value of scale, want <= 1, got %v", scale)
+
+ // Check any invalid modifications in read-only input.
+ if !equalGeneral(tl, tlCopy) {
+ t.Errorf("%v: unexpected modification in TL", name)
}
- if n1 == 0 || n2 == 0 {
- return
+ if !equalGeneral(tr, trCopy) {
+ t.Errorf("%v: unexpected modification in TR", name)
+ }
+ if !equalGeneral(b, bCopy) {
+ t.Errorf("%v: unexpected modification in B", name)
}
-
- prefix := fmt.Sprintf("Case n1=%v, n2=%v, isgn=%v", n1, n2, isgn)
// Check any invalid modifications of x.
if !generalOutsideAllNaN(x) {
- t.Errorf("%v: out-of-range write to x\n%v", prefix, x.Data)
+ t.Errorf("%v: out-of-range write to x\n%v", name, x.Data)
}
- var xnormWant float64
- for i := 0; i < n1; i++ {
- var rowsum float64
- for j := 0; j < n2; j++ {
- rowsum += math.Abs(x.Data[i*x.Stride+j])
- }
- if rowsum > xnormWant {
- xnormWant = rowsum
+ if n1 == 0 || n2 == 0 {
+ return
}
+
+ if scale <= 0 || 1 < scale {
+ t.Errorf("%v: invalid value of scale, want in (0,1], got %v", name, scale)
}
+
+ xnormWant := dlange(lapack.MaxRowSum, x.Rows, x.Cols, x.Data, x.Stride)
if xnormWant != xnorm {
- t.Errorf("%v: unexpected xnorm: want %v, got %v", prefix, xnormWant, xnorm)
+ t.Errorf("%v: unexpected xnorm: want %v, got %v", name, xnormWant, xnorm)
}
- // Multiply b by scale to get the wanted right-hand side.
- for i := 0; i < n1; i++ {
- for j := 0; j < n2; j++ {
- b.Data[i*b.Stride+j] *= scale
- }
+ if !ok {
+ t.Logf("%v: Dlasy2 returned ok=false", name)
+ return
}
- // Compute the wanted left-hand side.
- lhsWant := randomGeneral(n1, n2, n2, rnd)
+
+ // Compute diff := op(TL)*X + sgn*X*op(TR) - scale*B.
+ diff := zeros(n1, n2, n2)
if tranl {
- blas64.Gemm(blas.Trans, blas.NoTrans, 1, tl, x, 0, lhsWant)
+ blas64.Gemm(blas.Trans, blas.NoTrans, 1, tl, x, 0, diff)
} else {
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tl, x, 0, lhsWant)
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tl, x, 0, diff)
}
if tranr {
- blas64.Gemm(blas.NoTrans, blas.Trans, float64(isgn), x, tr, 1, lhsWant)
+ blas64.Gemm(blas.NoTrans, blas.Trans, float64(isgn), x, tr, 1, diff)
} else {
- blas64.Gemm(blas.NoTrans, blas.NoTrans, float64(isgn), x, tr, 1, lhsWant)
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, float64(isgn), x, tr, 1, diff)
}
- // Compare them.
for i := 0; i < n1; i++ {
for j := 0; j < n2; j++ {
- diff := lhsWant.Data[i*lhsWant.Stride+j] - b.Data[i*b.Stride+j]
- if math.Abs(diff) > tol && ok {
- t.Errorf("%v: unexpected result, diff[%v,%v]=%v", prefix, i, j, diff)
+ diff.Data[i*diff.Stride+j] -= scale * b.Data[i*b.Stride+j]
}
}
+ // Check that residual |op(TL)*X + sgn*X*op(TR) - scale*B| / |X| is small.
+ resid := dlange(lapack.MaxColumnSum, n1, n2, diff.Data, diff.Stride) / xnorm
+ if resid > tol {
+ t.Errorf("%v: unexpected result, resid=%v, want<=%v", name, resid, tol)
}
}

lapack/testlapack: extend test for Dlasy2 to cover rescaling