Repository: gonum/gonum
Pull request #1652 (2021-06-05 15:48:54)

lapack/gonum: add Dgesc2
So the checklist is NOT complete. I intend for this to be a first review of how LAPACK routines are implemented in Go and get some feedback. I'm sure my indexing is messed up, was quite confused seeeing how fortran alternated betweeen double and single indexing.

No tests implemented.
<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


soypat: Some lessons learned:
* Parenthesis matter. They can prevent precision loss and overflow if done correctly.
* `lda` is the "leading dimension" and used to index matrices, not `n`/`m` (actual dimensions)
* Documentation must be taken care of properly
* Add bounds/length checks not present in reference

As for the test, I need to generate not only the general matrix but also `ipiv` and `jpiv`. I'll take a look at other tests which do this i.e. `dlaqp2.go` and `dlaqps.go`

vladimir-ch: > As for the test, I need to generate not only the general matrix but also ipiv and jpiv.

It will be much easier (and more correct) if you first implement Dgetc2 and use it in the test to generate input for Dgesc2.
soypat: > It will be much easier (and more correct) if you first implement Dgetc2 and use it in the test to generate input for Dgesc2.

should I do it in this PR? should I start another?

vladimir-ch: `Dgetc2` should be a separate PR.
codecov-commenter: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
> Merging [#1652](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (504f670) into [master](https://codecov.io/gh/gonum/gonum/commit/e4c04dde5f994bd3aac25bb28fd5e796d9a728e4?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (e4c04dd) will **decrease** coverage by `0.18%`.
> The diff coverage is `39.39%`.

[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1652/graphs/tree.svg?width=650&height=150&src=pr&token=B7qHnqSPNL&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)

```diff
@@            Coverage Diff             @@
##           master    #1652      +/-   ##
==========================================
- Coverage   73.24%   73.05%   -0.19%     
==========================================
  Files         506      491      -15     
  Lines       59953    59463     -490     
==========================================
- Hits        43910    43442     -468     
+ Misses      13557    13551       -6     
+ Partials     2486     2470      -16     
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [lapack/gonum/dgesc2.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bGFwYWNrL2dvbnVtL2RnZXNjMi5nbw==) | `39.39% <39.39%> (ø)` | |
| [internal/asm/c64/conj.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2M2NC9jb25qLmdv) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [internal/asm/f64/scal.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9zY2FsLmdv) | | |
| [mat/offset\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L29mZnNldF9hcHBlbmdpbmUuZ28=) | | |
| [mat/index\_no\_bound\_checks.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L2luZGV4X25vX2JvdW5kX2NoZWNrcy5nbw==) | | |
| [internal/asm/f32/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2YzMi9zdHVic19ub2FzbS5nbw==) | | |
| [graph/iterator/nodes\_map\_safe.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaXRlcmF0b3Ivbm9kZXNfbWFwX3NhZmUuZ28=) | | |
| [...fourier/internal/fftpack/array\_no\_bounds\_checks.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-ZHNwL2ZvdXJpZXIvaW50ZXJuYWwvZmZ0cGFjay9hcnJheV9ub19ib3VuZHNfY2hlY2tzLmdv) | | |
| [internal/asm/c64/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2M2NC9zdHVic19ub2FzbS5nbw==) | | |
| [graph/internal/set/same\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaW50ZXJuYWwvc2V0L3NhbWVfYXBwZW5naW5lLmdv) | | |
| ... and [9 more](https://codecov.io/gh/gonum/gonum/pull/1652/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Last update [e4c04dd...504f670](https://codecov.io/gh/gonum/gonum/pull/1652?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).

soypat: Dgesc2 is passing albeit without the test you mentioned. I tried programming it and ended up leaving the `t.Error` declaration commented because it was failing although `X` was being calculated correctly.

Patch:
diff --git a/interp/cubic.go b/interp/cubic.go
index 4afb5e59..e6f0ed0a 100644
--- a/interp/cubic.go
+++ b/interp/cubic.go
@@ -377,8 +377,9 @@ func makeCubicSplineSecondDerivativeEquations(xs, ys []float64) (*mat.Tridiag, m
 
 // NaturalCubic is a piecewise cubic 1-dimensional interpolator with
 // continuous value, first and second derivatives, which can be fitted to (X, Y)
-// value pairs without providing derivatives. See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf
-// for details.
+// value pairs without providing derivatives. It uses the boundary conditions
+// Y''(left end ) = Y''(right end) = 0.
+// See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf for details.
 type NaturalCubic struct {
 	cubic PiecewiseCubic
 }
@@ -412,3 +413,48 @@ func (nc *NaturalCubic) Fit(xs, ys []float64) error {
 	}
 	return err
 }
+
+// ClampedCubic is a piecewise cubic 1-dimensional interpolator with
+// continuous value, first and second derivatives, which can be fitted to (X, Y)
+// value pairs without providing derivatives. It uses the boundary conditions
+// Y'(left end ) = Y'(right end) = 0.
+type ClampedCubic struct {
+	cubic PiecewiseCubic
+}
+
+// Predict returns the interpolation value at x.
+func (cc *ClampedCubic) Predict(x float64) float64 {
+	return cc.cubic.Predict(x)
+}
+
+// PredictDerivative returns the predicted derivative at x.
+func (cc *ClampedCubic) PredictDerivative(x float64) float64 {
+	return cc.cubic.PredictDerivative(x)
+}
+
+// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+// It panics if len(xs) < 2, elements of xs are not strictly increasing
+// or len(xs) != len(ys). It returns an error if solving the required system
+// of linear equations fails.
+func (cc *ClampedCubic) Fit(xs, ys []float64) error {
+	a, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+	// Add boundary conditions y''(left) = y''(right) = 0:
+	n := len(xs)
+	// Condition Y'(left end) = 0:
+	dxL := xs[1] - xs[0]
+	b.SetVec(0, (ys[1]-ys[0])/dxL)
+	a.SetBand(0, 0, dxL/3)
+	a.SetBand(0, 1, dxL/6)
+	// Condition Y'(right end) = 0:
+	m := n - 1
+	dxR := xs[m] - xs[m-1]
+	b.SetVec(m, (ys[m]-ys[m-1])/dxR)
+	a.SetBand(m, m, -dxR/3)
+	a.SetBand(m, m-1, -dxR/6)
+	x := mat.NewVecDense(n, nil)
+	err := a.SolveVecTo(x, false, b)
+	if err == nil {
+		cc.cubic.fitWithSecondDerivatives(xs, ys, x.RawVector().Data)
+	}
+	return err
+}
diff --git a/lapack/gonum/dgesc2.go b/lapack/gonum/dgesc2.go
new file mode 100644
index 00000000..4d0f1260
--- /dev/null
+++ b/lapack/gonum/dgesc2.go
@@ -0,0 +1,83 @@
+// Copyright ©2021 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 gonum
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/blas/blas64"
+)
+
+// Dgesc2 solves a system of linear equations
+//  A * X = scale * RHS
+// with a general N-by-N matrix A using the LU factorization with
+// complete pivoting computed by Dgetc2. The result is placed in
+// rhs on exit.
+//
+// Dgesc2 is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dgesc2(n int, a []float64, lda int, rhs []float64, ipiv, jpiv []int) (scale float64) {
+	switch {
+	case n < 0:
+		panic(nLT0)
+	case lda < max(1, n):
+		panic(badLdA)
+	}
+
+	// Quick return if possible.
+	if n == 0 {
+		return 0
+	}
+
+	switch {
+	case len(a) < (n-1)*lda+n:
+		panic(shortA)
+	case len(rhs) < n:
+		panic(shortRHS)
+	case len(ipiv) != n:
+		panic(badLenIpiv)
+	case len(jpiv) != n:
+		panic(badLenJpiv)
+	}
+
+	const smlnum = dlamchS / dlamchP
+	if len(a) < (n-1)*lda+n {
+		panic(shortA)
+	}
+
+	// Apply permutations ipiv to RHS.
+	impl.Dlaswp(1, rhs, 1, 0, n-1, ipiv[:n], 1)
+
+	// Solve for L part.
+	for i := 0; i < n-1; i++ {
+		for j := i + 1; j < n; j++ {
+			rhs[j] -= float64(a[j*lda+i] * rhs[i])
+		}
+	}
+
+	// Solve for U part.
+
+	scale = 1.0
+
+	// Check for scaling.
+	bi := blas64.Implementation()
+	i := bi.Idamax(n, rhs, 1)
+	if 2*smlnum*math.Abs(rhs[i]) > math.Abs(a[(n-1)*lda+(n-1)]) {
+		temp := 0.5 / math.Abs(rhs[i])
+		bi.Dscal(n, temp, rhs, 1)
+		scale *= temp
+	}
+
+	for i := n - 1; i >= 0; i-- {
+		temp := 1.0 / a[i*lda+i]
+		rhs[i] *= temp
+		for j := i + 1; j < n; j++ {
+			rhs[i] -= float64(rhs[j] * (a[i*lda+j] * temp))
+		}
+	}
+
+	// Apply permutations jpiv to the solution (rhs).
+	impl.Dlaswp(1, rhs, 1, 0, n-1, jpiv[:n], -1)
+	return scale
+}
diff --git a/lapack/gonum/dormbr.go b/lapack/gonum/dormbr.go
index e1b37373..20f3ca57 100644
--- a/lapack/gonum/dormbr.go
+++ b/lapack/gonum/dormbr.go
@@ -145,7 +145,7 @@ func (impl Implementation) Dormbr(vect lapack.ApplyOrtho, side blas.Side, trans
 				i1 = 1
 				i2 = 0
 			}
-			impl.Dormqr(side, trans, mi, ni, nq-1, a[1*lda:], lda, tau[:nq-1], c[i1*ldc+i2:], ldc, work, lwork)
+			impl.Dormqr(side, trans, mi, ni, nq-1, a[lda:], lda, tau[:nq-1], c[i1*ldc+i2:], ldc, work, lwork)
 		}
 		work[0] = float64(lworkopt)
 		return
diff --git a/lapack/gonum/errors.go b/lapack/gonum/errors.go
index 28f3b150..1429e65d 100644
--- a/lapack/gonum/errors.go
+++ b/lapack/gonum/errors.go
@@ -101,6 +101,7 @@ const (
 	badLenAlpha    = "lapack: bad length of alpha"
 	badLenBeta     = "lapack: bad length of beta"
 	badLenIpiv     = "lapack: bad length of ipiv"
+	badLenJpiv     = "lapack: bad length of jpiv"
 	badLenJpvt     = "lapack: bad length of jpvt"
 	badLenK        = "lapack: bad length of k"
 	badLenSelected = "lapack: bad length of selected"
@@ -126,6 +127,7 @@ const (
 	shortIWork = "lapack: insufficient length of iwork"
 	shortIsgn  = "lapack: insufficient length of isgn"
 	shortQ     = "lapack: insufficient length of q"
+	shortRHS   = "lapack: insufficient length of rhs"
 	shortS     = "lapack: insufficient length of s"
 	shortScale = "lapack: insufficient length of scale"
 	shortT     = "lapack: insufficient length of t"


Test patch:
diff --git a/interp/cubic_test.go b/interp/cubic_test.go
index b59df3e8..63b24cb1 100644
--- a/interp/cubic_test.go
+++ b/interp/cubic_test.go
@@ -799,6 +799,48 @@ func TestNaturalCubicFit(t *testing.T) {
 		if math.Abs(got-want[i]) > 1e-14 {
 			t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
 		}
+		got = nc.PredictDerivative(testXs[i])
+		want := discrDerivPredict(&nc, xs[0], xs[len(xs)-1], testXs[i], 1e-9)
+		if math.Abs(got-want) > 1e-6 {
+			t.Errorf("Mismatch in predicted dY/dX value for x = %g: got %v, want %g", testXs[i], got, want)
+		}
 	}
+}
 
+func TestClampedCubicFit(t *testing.T) {
+	t.Parallel()
+	xs := []float64{-1, 0, 2, 3.5}
+	ys := []float64{2, 0, 2, 1.5}
+	var cc ClampedCubic
+	err := cc.Fit(xs, ys)
+	if err != nil {
+		t.Errorf("Error when fitting ClampedCubic: %v", err)
+	}
+	testXs := []float64{-1, -0.99, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.49, 3.5}
+	// From scipy.interpolate.CubicSpline:
+	want := []float64{
+		2.0,
+		1.999564111111111,
+		1.2021604938271606,
+		0.0,
+		-0.20833333333333343,
+		0.41975308641975284,
+		1.337962962962962,
+		2.0,
+		2.026748971193416,
+		1.7078189300411522,
+		1.5001129711934151,
+		1.4999999999999996,
+	}
+	for i := 0; i < len(testXs); i++ {
+		got := cc.Predict(testXs[i])
+		if math.Abs(got-want[i]) > 1e-14 {
+			t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
+		}
+		got = cc.PredictDerivative(testXs[i])
+		want := discrDerivPredict(&cc, xs[0], xs[len(xs)-1], testXs[i], 1e-9)
+		if math.Abs(got-want) > 1e-6 {
+			t.Errorf("Mismatch in predicted dY/dX value for x = %g: got %v, want %g", testXs[i], got, want)
+		}
+	}
 }
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index 2dded735..a296d64f 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -92,6 +92,11 @@ func TestDgerq2(t *testing.T) {
 	testlapack.Dgerq2Test(t, impl)
 }
 
+func TestDgesc2(t *testing.T) {
+	t.Parallel()
+	testlapack.Dgesc2Test(t, impl)
+}
+
 func TestDgeqp3(t *testing.T) {
 	t.Parallel()
 	testlapack.Dgeqp3Test(t, impl)
diff --git a/lapack/testlapack/dgesc2.go b/lapack/testlapack/dgesc2.go
new file mode 100644
index 00000000..4c351b51
--- /dev/null
+++ b/lapack/testlapack/dgesc2.go
@@ -0,0 +1,91 @@
+// Copyright ©2021 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 Dgesc2er interface {
+	Dgetc2er
+	// Dgesc2 solves a system of linear equations
+	//  A * X = scale * RHS
+	// with a general n×n matrix A using the LU factorization with
+	// complete pivoting computed by Dgetc2. The result is placed in
+	// rhs on exit.
+	Dgesc2(n int, a []float64, lda int, rhs []float64, ipiv, jpiv []int) (scale float64)
+}
+
+func Dgesc2Test(t *testing.T, impl Dgesc2er) {
+	const tol = 1e-12
+	rnd := rand.New(rand.NewSource(1))
+	for _, test := range []struct {
+		n, lda int
+	}{
+		{3, 0},
+		{5, 0},
+		{20, 30},
+		{200, 0},
+	} {
+		testSolveDgesc2(t, impl, rnd, test.n, test.lda, tol)
+	}
+}
+
+func testSolveDgesc2(t *testing.T, impl Dgesc2er, rnd *rand.Rand, n, lda int, tol float64) {
+	name := fmt.Sprintf("n=%v,lda=%v", n, lda)
+	if lda == 0 {
+		lda = n
+	}
+	// Generate random general matrix.
+	a := randomGeneral(n, n, lda, rnd)
+	// anorm := floats.Norm(a.Data, 1)
+
+	// Generate a random solution.
+	xWant := randomGeneral(n, 1, 1, rnd)
+	// xnorm := floats.Norm(xWant.Data, 1)
+
+	// Compute RHS vector that solves for X such that  A*X = scale * RHS
+	rhs := zeros(n, 1, 1)
+	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, a, xWant, 1, rhs)
+	rhsCopy := zeros(n, 1, 1) // Will contain A*x result.
+	copyGeneral(rhsCopy, rhs)
+	// Compute LU factorization with full pivoting.
+	lu := zeros(n, n, lda)
+	copyGeneral(lu, a)
+	ipiv := make([]int, n)
+	jpiv := make([]int, n)
+	impl.Dgetc2(n, lu.Data, lu.Stride, ipiv, jpiv)
+
+	// Solve using lu factorization.
+	scale := impl.Dgesc2(lu.Rows, lu.Data, lu.Stride, rhs.Data, ipiv, jpiv)
+	x := rhs
+	if scale < 0 || scale > 1 {
+		t.Errorf("%v: resulting scale out of bounds [0,1]", name)
+	}
+
+	var diff float64
+	for i := range x.Data {
+		diff = math.Max(diff, math.Abs(xWant.Data[i]-x.Data[i]))
+	}
+	if diff > tol {
+		t.Errorf("%v: unexpected result, diff=%v", name, diff)
+	}
+	// |A*X - scale*RHS| / |A| / |X| is an indicator that solution is good
+	// AxResult := zeros(n, 1, 1)
+	// blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, a, x, 1, AxResult)
+	// blas64.Scal(scale, blas64.Vector{N: n, Data: rhsCopy.Data, Inc: 1})
+	// floats.Sub(AxResult.Data, rhsCopy.Data)
+
+	// residualNorm := floats.Norm(rhsCopy.Data, 1) / anorm / xnorm
+	// if residualNorm > tol {
+	// 	t.Errorf("%v: |A*X - scale*RHS| / |A| / |X| = %g is greater than permissible tol", name, residualNorm)
+	// }
+}
diff --git a/lapack/testlapack/dlasy2.go b/lapack/testlapack/dlasy2.go
index 1128fd95..6e1a1f6f 100644
--- a/lapack/testlapack/dlasy2.go
+++ b/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,9 +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)
 							}
 						}
 					}
@@ -38,68 +42,83 @@ func Dlasy2Test(t *testing.T, impl Dlasy2er) {
 	}
 }
 
-func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, rnd *rand.Rand) {
-	const tol = 1e-10
+func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, big bool, rnd *rand.Rand) {
+	const tol = 1e-14
+
+	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)
+	}
 }