Repository: gonum/gonum
Pull request #1603 (2021-04-09 10:39:49)

lapack/gonum: check for overflow in Dtgsja
Please take a look.

I don't have a failing test case for this change but it was done in the reference in https://github.com/Reference-LAPACK/lapack/pull/502 and looks like a safe change to make. The PR references a test case for CGGSVD3 which does trigger a failure but it looks like a deficiency of our test (too strict by not taking into account the norm of the input matrices).

<!--
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 _____."
-->


kortschak: it will be nice when fuzzing is built in and there is support for structured data. The first part of that is coming and the second part is proposed. With that, finding test cases will be easier for things like this.

Patch:
diff --git a/lapack/gonum/dtgsja.go b/lapack/gonum/dtgsja.go
index 3628ca62..5bc0ac1a 100644
--- a/lapack/gonum/dtgsja.go
+++ b/lapack/gonum/dtgsja.go
@@ -326,10 +326,8 @@ func (impl Implementation) Dtgsja(jobU, jobV, jobQ lapack.GSVDJob, m, p, n, k, l
 				for i := 0; i < min(l, m-k); i++ {
 					a1 := a[(k+i)*lda+n-l+i]
 					b1 := b[i*ldb+n-l+i]
-
-					if a1 != 0 {
-						gamma := b1 / a1
-
+					gamma := b1 / a1
+					if !math.IsInf(gamma, 0) {
 						// Change sign if necessary.
 						if gamma < 0 {
 							bi.Dscal(l-i, -1, b[i*ldb+n-l+i:], 1)
diff --git a/mat/band.go b/mat/band.go
index 5654f906..f2fbed06 100644
--- a/mat/band.go
+++ b/mat/band.go
@@ -134,10 +134,10 @@ func NewBandDense(r, c, kl, ku int, data []float64) *BandDense {
 		if r == 0 || c == 0 {
 			panic(ErrZeroLength)
 		}
-		panic("mat: negative dimension")
+		panic(ErrNegativeDimension)
 	}
 	if kl+1 > r || ku+1 > c {
-		panic("mat: band out of range")
+		panic(ErrBandwidth)
 	}
 	bc := kl + ku + 1
 	if data != nil && len(data) != min(r, c+kl)*bc {
diff --git a/mat/index_bound_checks.go b/mat/index_bound_checks.go
index 59815a67..cb4c60bd 100644
--- a/mat/index_bound_checks.go
+++ b/mat/index_bound_checks.go
@@ -346,3 +346,52 @@ func (d *DiagDense) setDiag(i int, v float64) {
 	}
 	d.mat.Data[i*d.mat.Inc] = v
 }
+
+// At returns the element at row i, column j.
+func (a *Tridiag) At(i, j int) float64 {
+	return a.at(i, j)
+}
+
+func (a *Tridiag) at(i, j int) float64 {
+	if uint(i) >= uint(d.mat.N) {
+		panic(ErrRowAccess)
+	}
+	if uint(j) >= uint(d.mat.N) {
+		panic(ErrColAccess)
+	}
+	switch i - j {
+	case -1:
+		return a.mat.DU[i]
+	case 0:
+		return a.mat.D[i]
+	case 1:
+		return a.mat.DL[j]
+	default:
+		return 0
+	}
+}
+
+// SetBand sets the element at row i, column j to the value v.
+// It panics if the location is outside the appropriate region of the matrix.
+func (a *Tridiag) SetBand(i, j int, v float64) {
+	a.set(i, j, v)
+}
+
+func (a *Tridiag) set(i, j int, v float64) {
+	if uint(i) >= uint(a.mat.N) {
+		panic(ErrRowAccess)
+	}
+	if uint(j) >= uint(a.mat.N) {
+		panic(ErrColAccess)
+	}
+	switch i - j {
+	case -1:
+		a.mat.DU[i] = v
+	case 0:
+		a.mat.D[i] = v
+	case 1:
+		a.mat.DL[j] = v
+	default:
+		panic(ErrBandSet)
+	}
+}
diff --git a/mat/index_no_bound_checks.go b/mat/index_no_bound_checks.go
index 051f8437..a15228a0 100644
--- a/mat/index_no_bound_checks.go
+++ b/mat/index_no_bound_checks.go
@@ -357,3 +357,52 @@ func (d *DiagDense) SetDiag(i int, v float64) {
 func (d *DiagDense) setDiag(i int, v float64) {
 	d.mat.Data[i*d.mat.Inc] = v
 }
+
+// At returns the element at row i, column j.
+func (a *Tridiag) At(i, j int) float64 {
+	if uint(i) >= uint(a.mat.N) {
+		panic(ErrRowAccess)
+	}
+	if uint(j) >= uint(a.mat.N) {
+		panic(ErrColAccess)
+	}
+	return a.at(i, j)
+}
+
+func (a *Tridiag) at(i, j int) float64 {
+	switch i - j {
+	case -1:
+		return a.mat.DU[i]
+	case 0:
+		return a.mat.D[i]
+	case 1:
+		return a.mat.DL[j]
+	default:
+		return 0
+	}
+}
+
+// SetBand sets the element at row i, column j to the value v.
+// It panics if the location is outside the appropriate region of the matrix.
+func (a *Tridiag) SetBand(i, j int, v float64) {
+	if uint(i) >= uint(a.mat.N) {
+		panic(ErrRowAccess)
+	}
+	if uint(j) >= uint(a.mat.N) {
+		panic(ErrColAccess)
+	}
+	a.set(i, j, v)
+}
+
+func (a *Tridiag) set(i, j int, v float64) {
+	switch i - j {
+	case -1:
+		a.mat.DU[i] = v
+	case 0:
+		a.mat.D[i] = v
+	case 1:
+		a.mat.DL[j] = v
+	default:
+		panic(ErrBandSet)
+	}
+}
diff --git a/mat/matrix.go b/mat/matrix.go
index 799a173a..98c34612 100644
--- a/mat/matrix.go
+++ b/mat/matrix.go
@@ -244,7 +244,7 @@ func untranspose(a Matrix) (Matrix, bool) {
 func untransposeExtract(a Matrix) (Matrix, bool) {
 	ut, trans := untranspose(a)
 	switch m := ut.(type) {
-	case *DiagDense, *SymBandDense, *TriBandDense, *BandDense, *TriDense, *SymDense, *Dense, *VecDense:
+	case *DiagDense, *SymBandDense, *TriBandDense, *BandDense, *TriDense, *SymDense, *Dense, *VecDense, *Tridiag:
 		return m, trans
 	// TODO(btracey): Add here if we ever have an equivalent of RawDiagDense.
 	case RawSymBander:
@@ -291,6 +291,10 @@ func untransposeExtract(a Matrix) (Matrix, bool) {
 		var v VecDense
 		v.SetRawVector(m.RawVector())
 		return &v, trans
+	case RawTridiagonaler:
+		var d Tridiag
+		d.SetRawTridiagonal(m.RawTridiagonal())
+		return &d, trans
 	default:
 		return ut, trans
 	}
diff --git a/mat/tridiag.go b/mat/tridiag.go
new file mode 100644
index 00000000..afac4846
--- /dev/null
+++ b/mat/tridiag.go
@@ -0,0 +1,397 @@
+// Copyright ©2020 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 mat
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/blas"
+	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/internal/asm/f64"
+	"gonum.org/v1/gonum/lapack/lapack64"
+)
+
+var (
+	tridiagDense *Tridiag
+	_            Matrix           = tridiagDense
+	_            allMatrix        = tridiagDense
+	_            denseMatrix      = tridiagDense
+	_            Banded           = tridiagDense
+	_            MutableBanded    = tridiagDense
+	_            RawTridiagonaler = tridiagDense
+)
+
+// A RawTridiagonaler can return a lapack64.Tridiagonal representation of the
+// receiver. Changes to the elements of DL, D, DU in lapack64.Tridiagonal will
+// be reflected in the original matrix, changes to the N field will not.
+type RawTridiagonaler interface {
+	RawTridiagonal() lapack64.Tridiagonal
+}
+
+// Tridiag represents a tridiagonal matrix by its three diagonals.
+type Tridiag struct {
+	mat lapack64.Tridiagonal
+}
+
+// NewTridiag creates a new n×n tridiagonal matrix with the first sub-diagonal
+// in dl, the main diagonal in d and the first super-diagonal in du. If all of
+// dl, d, and du are nil, new backing slices will be allocated for them. If dl
+// and du have length n-1 and d has length n, they will be used as backing
+// slices, and changes to the elements of the returned Tridiag will be reflected
+// in dl, d, du. If neither of these is true, NewTridiag will panic.
+func NewTridiag(n int, dl, d, du []float64) *Tridiag {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
+		panic(ErrNegativeDimension)
+	}
+	if dl != nil || d != nil || du != nil {
+		if len(dl) != n-1 || len(d) != n || len(du) != n-1 {
+			panic(ErrShape)
+		}
+	} else {
+		d = make([]float64, n)
+		if n > 1 {
+			dl = make([]float64, n-1)
+			du = make([]float64, n-1)
+		}
+	}
+	return &Tridiag{
+		mat: lapack64.Tridiagonal{
+			N:  n,
+			DL: dl,
+			D:  d,
+			DU: du,
+		},
+	}
+}
+
+// Dims returns the number of rows and columns in the matrix.
+func (a *Tridiag) Dims() (r, c int) {
+	return a.mat.N, a.mat.N
+}
+
+// Bandwidth returns 1, 1 - the upper and lower bandwidths of the matrix.
+func (a *Tridiag) Bandwidth() (kl, ku int) {
+	return 1, 1
+}
+
+// T performs an implicit transpose by returning the receiver inside a Transpose.
+func (a *Tridiag) T() Matrix {
+	// An alternative would be to return the receiver with DL,DU swapped; the
+	// untranspose function would then always return false. With Transpose the
+	// diagonal swapping will be done in tridiagonal routines in lapack like
+	// lapack64.Gtsv or gonum.Dlagtm based on the trans parameter.
+	return Transpose{a}
+}
+
+// TBand performs an implicit transpose by returning the receiver inside a
+// TransposeBand.
+func (a *Tridiag) TBand() Banded {
+	// An alternative would be to return the receiver with DL,DU swapped; see
+	// explanation in T above.
+	return TransposeBand{a}
+}
+
+// RawTridiagonal returns the underlying lapack64.Tridiagonal used by the
+// receiver. Changes to elements in the receiver following the call will be
+// reflected in the returned matrix.
+func (a *Tridiag) RawTridiagonal() lapack64.Tridiagonal {
+	return a.mat
+}
+
+// SetRawTridiagonal sets the underlying lapack64.Tridiagonal used by the
+// receiver. Changes to elements in the receiver following the call will be
+// reflected in the input.
+func (a *Tridiag) SetRawTridiagonal(mat lapack64.Tridiagonal) {
+	a.mat = mat
+}
+
+// IsEmpty returns whether the receiver is empty. Empty matrices can be the
+// receiver for size-restricted operations. The receiver can be zeroed using
+// Reset.
+func (a *Tridiag) IsEmpty() bool {
+	return a.mat.N == 0
+}
+
+// Reset empties the matrix so that it can be reused as the receiver of a
+// dimensionally restricted operation.
+//
+// Reset should not be used when the matrix shares backing data. See the Reseter
+// interface for more information.
+func (a *Tridiag) Reset() {
+	a.mat.N = 0
+	a.mat.DL = a.mat.DL[:0]
+	a.mat.D = a.mat.D[:0]
+	a.mat.DU = a.mat.DU[:0]
+}
+
+// Clone makes a copy of the input Tridiag into the receiver, overwriting the
+// previous value of the receiver. Clone does not place any restrictions on
+// receiver shape.
+func (a *Tridiag) Clone(from *Tridiag) {
+	n := from.mat.N
+	switch n {
+	case 0:
+		panic(ErrZeroLength)
+	case 1:
+		a.mat = lapack64.Tridiagonal{
+			N:  1,
+			DL: use(a.mat.DL, 0),
+			D:  use(a.mat.D, 1),
+			DU: use(a.mat.DU, 0),
+		}
+		a.mat.D[0] = from.mat.D[0]
+	default:
+		a.mat = lapack64.Tridiagonal{
+			N:  n,
+			DL: use(a.mat.DL, n-1),
+			D:  use(a.mat.D, n),
+			DU: use(a.mat.DU, n-1),
+		}
+		copy(a.mat.DL, from.mat.DL)
+		copy(a.mat.D, from.mat.D)
+		copy(a.mat.DU, from.mat.DU)
+	}
+}
+
+// DiagView returns the diagonal as a matrix backed by the original data.
+func (a *Tridiag) DiagView() Diagonal {
+	return &DiagDense{
+		mat: blas64.Vector{
+			N:    a.mat.N,
+			Data: a.mat.D[:a.mat.N],
+			Inc:  1,
+		},
+	}
+}
+
+// Zero sets all of the matrix elements to zero.
+func (a *Tridiag) Zero() {
+	zero(a.mat.DL)
+	zero(a.mat.D)
+	zero(a.mat.DU)
+}
+
+// Trace returns the trace of the matrix.
+func (a *Tridiag) Trace() float64 {
+	return f64.Sum(a.mat.D)
+}
+
+// Norm returns the specified norm of the receiver. Valid norms are:
+//  1 - The maximum absolute column sum
+//  2 - The Frobenius norm, the square root of the sum of the squares of the elements
+//  Inf - The maximum absolute row sum
+//
+// Norm will panic with ErrNormOrder if an illegal norm is specified.
+func (a *Tridiag) Norm(norm float64) float64 {
+	return lapack64.Langt(normLapack(norm, false), a.mat)
+}
+
+// MulVecTo computes A⋅x or Aᵀ⋅x storing the result into dst.
+func (a *Tridiag) MulVecTo(dst *VecDense, trans bool, x Vector) {
+	n := a.mat.N
+	if x.Len() != n {
+		panic(ErrShape)
+	}
+	dst.reuseAsNonZeroed(n)
+	t := blas.NoTrans
+	if trans {
+		t = blas.Trans
+	}
+	xMat, _ := untransposeExtract(x)
+	if xVec, ok := xMat.(*VecDense); ok && dst != xVec {
+		dst.checkOverlap(xVec.mat)
+		lapack64.Lagtm(t, 1, a.mat, xVec.asGeneral(), 0, dst.asGeneral())
+	} else {
+		xCopy := getWorkspaceVec(n, false)
+		xCopy.CloneFromVec(x)
+		lapack64.Lagtm(t, 1, a.mat, xCopy.asGeneral(), 0, dst.asGeneral())
+		putWorkspaceVec(xCopy)
+	}
+}
+
+// SolveTo solves a tridiagonal system A⋅X = B  or  Aᵀ⋅X = B where A is an
+// n×n tridiagonal matrix represented by the receiver and B is a given n×nrhs
+// matrix. If A is non-singular, the result will be stored into dst and nil will
+// be returned. If A is singular, the contents of dst will be undefined and a
+// Condition error will be returned.
+func (a *Tridiag) SolveTo(dst *Dense, trans bool, b Matrix) error {
+	n, nrhs := b.Dims()
+	if n != a.mat.N {
+		panic(ErrShape)
+	}
+	if b, ok := b.(RawMatrixer); ok && dst != b {
+		dst.checkOverlap(b.RawMatrix())
+	}
+	dst.reuseAsNonZeroed(n, nrhs)
+	if dst != b {
+		dst.Copy(b)
+	}
+	var aCopy Tridiag
+	aCopy.Clone(a)
+	var ok bool
+	if trans {
+		ok = lapack64.Gtsv(blas.Trans, aCopy.mat, dst.mat)
+	} else {
+		ok = lapack64.Gtsv(blas.NoTrans, aCopy.mat, dst.mat)
+	}
+	if !ok {
+		return Condition(math.Inf(1))
+	}
+	return nil
+}
+
+// SolveVecTo solves a tridiagonal system A⋅X = B  or  Aᵀ⋅X = B where A is an
+// n×n tridiagonal matrix represented by the receiver and b is a given n-vector.
+// If A is non-singular, the result will be stored into dst and nil will be
+// returned. If A is singular, the contents of dst will be undefined and a
+// Condition error will be returned.
+func (a *Tridiag) SolveVecTo(dst *VecDense, trans bool, b Vector) error {
+	n, nrhs := b.Dims()
+	if n != a.mat.N || nrhs != 1 {
+		panic(ErrShape)
+	}
+	if b, ok := b.(RawVectorer); ok && dst != b {
+		dst.checkOverlap(b.RawVector())
+	}
+	dst.reuseAsNonZeroed(n)
+	if dst != b {
+		dst.CopyVec(b)
+	}
+	var aCopy Tridiag
+	aCopy.Clone(a)
+	var ok bool
+	if trans {
+		ok = lapack64.Gtsv(blas.Trans, aCopy.mat, dst.asGeneral())
+	} else {
+		ok = lapack64.Gtsv(blas.NoTrans, aCopy.mat, dst.asGeneral())
+	}
+	if !ok {
+		return Condition(math.Inf(1))
+	}
+	return nil
+}
+
+// DoNonZero calls the function fn for each of the non-zero elements of A. The
+// function fn takes a row/column index and the element value of A at (i,j).
+func (a *Tridiag) DoNonZero(fn func(i, j int, v float64)) {
+	for i, aij := range a.mat.DU {
+		if aij != 0 {
+			fn(i, i+1, aij)
+		}
+	}
+	for i, aii := range a.mat.D {
+		if aii != 0 {
+			fn(i, i, aii)
+		}
+	}
+	for i, aij := range a.mat.DL {
+		if aij != 0 {
+			fn(i+1, i, aij)
+		}
+	}
+}
+
+// DoRowNonZero calls the function fn for each of the non-zero elements of row i
+// of A. The function fn takes a row/column index and the element value of A at
+// (i,j).
+func (a *Tridiag) DoRowNonZero(i int, fn func(i, j int, v float64)) {
+	n := a.mat.N
+	if uint(i) >= uint(n) {
+		panic(ErrRowAccess)
+	}
+	if n == 1 {
+		v := a.mat.D[0]
+		if v != 0 {
+			fn(0, 0, v)
+		}
+		return
+	}
+	switch i {
+	case 0:
+		v := a.mat.D[0]
+		if v != 0 {
+			fn(i, 0, v)
+		}
+		v = a.mat.DU[0]
+		if v != 0 {
+			fn(i, 1, v)
+		}
+	case n - 1:
+		v := a.mat.DL[n-2]
+		if v != 0 {
+			fn(n-1, n-2, v)
+		}
+		v = a.mat.D[n-1]
+		if v != 0 {
+			fn(n-1, n-1, v)
+		}
+	default:
+		v := a.mat.DL[i-1]
+		if v != 0 {
+			fn(i, i-1, v)
+		}
+		v = a.mat.D[i]
+		if v != 0 {
+			fn(i, i, v)
+		}
+		v = a.mat.DU[i]
+		if v != 0 {
+			fn(i, i+1, v)
+		}
+	}
+}
+
+// DoColNonZero calls the function fn for each of the non-zero elements of
+// column j of A. The function fn takes a row/column index and the element value
+// of A at (i, j).
+func (a *Tridiag) DoColNonZero(j int, fn func(i, j int, v float64)) {
+	n := a.mat.N
+	if uint(j) >= uint(n) {
+		panic(ErrColAccess)
+	}
+	if n == 1 {
+		v := a.mat.D[0]
+		if v != 0 {
+			fn(0, 0, v)
+		}
+		return
+	}
+	switch j {
+	case 0:
+		v := a.mat.D[0]
+		if v != 0 {
+			fn(0, 0, v)
+		}
+		v = a.mat.DL[0]
+		if v != 0 {
+			fn(1, 0, v)
+		}
+	case n - 1:
+		v := a.mat.DU[n-2]
+		if v != 0 {
+			fn(n-2, n-1, v)
+		}
+		v = a.mat.D[n-1]
+		if v != 0 {
+			fn(n-1, n-1, v)
+		}
+	default:
+		v := a.mat.DU[j-1]
+		if v != 0 {
+			fn(j-1, j, v)
+		}
+		v = a.mat.D[j]
+		if v != 0 {
+			fn(j, j, v)
+		}
+		v = a.mat.DL[j]
+		if v != 0 {
+			fn(j+1, j, v)
+		}
+	}
+}


Test patch:
diff --git a/mat/basictypes_test.go b/mat/basictypes_test.go
index a881f0f2..22970e6a 100644
--- a/mat/basictypes_test.go
+++ b/mat/basictypes_test.go
@@ -23,6 +23,14 @@ func (m *basicMatrix) At(r, c int) float64 { return (*Dense)(m).At(r, c) }
 func (m *basicMatrix) Dims() (r, c int)    { return (*Dense)(m).Dims() }
 func (m *basicMatrix) T() Matrix           { return Transpose{m} }
 
+type rawMatrix struct {
+	*basicMatrix
+}
+
+func (a *rawMatrix) RawMatrix() blas64.General {
+	return a.basicMatrix.mat
+}
+
 type basicVector VecDense
 
 var _ Vector = &basicVector{}
diff --git a/mat/matrix_test.go b/mat/matrix_test.go
index 3d0dee09..b9131fb6 100644
--- a/mat/matrix_test.go
+++ b/mat/matrix_test.go
@@ -612,6 +612,12 @@ func TestDoer(t *testing.T) {
 		NewSymBandDense(3, 1, ones(3*(1+1))),
 		NewSymBandDense(6, 1, ones(6*(1+1))),
 		NewSymBandDense(6, 2, ones(6*(2+1))),
+		NewTridiag(1, nil, ones(1), nil),
+		NewTridiag(2, ones(1), ones(2), ones(1)),
+		NewTridiag(3, ones(2), ones(3), ones(2)),
+		NewTridiag(4, ones(3), ones(4), ones(3)),
+		NewTridiag(7, ones(6), ones(7), ones(6)),
+		NewTridiag(10, ones(9), ones(10), ones(9)),
 	} {
 		r, c := m.Dims()
 
@@ -714,6 +720,12 @@ func TestMulVecToer(t *testing.T) {
 		NewSymBandDense(10, 0, random(10)),
 		NewSymBandDense(10, 1, random(20)),
 		NewSymBandDense(10, 4, random(50)),
+		NewTridiag(1, nil, random(1), nil),
+		NewTridiag(2, random(1), random(2), random(1)),
+		NewTridiag(3, random(2), random(3), random(2)),
+		NewTridiag(4, random(3), random(4), random(3)),
+		NewTridiag(7, random(6), random(7), random(6)),
+		NewTridiag(10, random(9), random(10), random(9)),
 	} {
 		// Dense copy of A used for computing the expected result.
 		var aDense Dense
diff --git a/mat/tridiag_test.go b/mat/tridiag_test.go
new file mode 100644
index 00000000..2b385704
--- /dev/null
+++ b/mat/tridiag_test.go
@@ -0,0 +1,530 @@
+// 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 mat
+
+import (
+	"fmt"
+	"math"
+	"reflect"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/lapack/lapack64"
+)
+
+func TestNewTridiag(t *testing.T) {
+	for i, test := range []struct {
+		n         int
+		dl, d, du []float64
+		panics    bool
+		want      *Tridiag
+		dense     *Dense
+	}{
+		{
+			n:      1,
+			dl:     nil,
+			d:      []float64{1.2},
+			du:     nil,
+			panics: false,
+			want: &Tridiag{
+				mat: lapack64.Tridiagonal{
+					N:  1,
+					DL: nil,
+					D:  []float64{1.2},
+					DU: nil,
+				},
+			},
+			dense: NewDense(1, 1, []float64{1.2}),
+		},
+		{
+			n:      1,
+			dl:     []float64{},
+			d:      []float64{1.2},
+			du:     []float64{},
+			panics: false,
+			want: &Tridiag{
+				mat: lapack64.Tridiagonal{
+					N:  1,
+					DL: []float64{},
+					D:  []float64{1.2},
+					DU: []float64{},
+				},
+			},
+			dense: NewDense(1, 1, []float64{1.2}),
+		},
+		{
+			n:      4,
+			dl:     []float64{1.2, 2.3, 3.4},
+			d:      []float64{4.5, 5.6, 6.7, 7.8},
+			du:     []float64{8.9, 9.0, 0.1},
+			panics: false,
+			want: &Tridiag{
+				mat: lapack64.Tridiagonal{
+					N:  4,
+					DL: []float64{1.2, 2.3, 3.4},
+					D:  []float64{4.5, 5.6, 6.7, 7.8},
+					DU: []float64{8.9, 9.0, 0.1},
+				},
+			},
+			dense: NewDense(4, 4, []float64{
+				4.5, 8.9, 0, 0,
+				1.2, 5.6, 9.0, 0,
+				0, 2.3, 6.7, 0.1,
+				0, 0, 3.4, 7.8,
+			}),
+		},
+		{
+			n:      4,
+			dl:     nil,
+			d:      nil,
+			du:     nil,
+			panics: false,
+			want: &Tridiag{
+				mat: lapack64.Tridiagonal{
+					N:  4,
+					DL: []float64{0, 0, 0},
+					D:  []float64{0, 0, 0, 0},
+					DU: []float64{0, 0, 0},
+				},
+			},
+			dense: NewDense(4, 4, nil),
+		},
+		{
+			n:      -1,
+			panics: true,
+		},
+		{
+			n:      0,
+			panics: true,
+		},
+		{
+			n:      1,
+			dl:     []float64{1.2},
+			d:      nil,
+			du:     nil,
+			panics: true,
+		},
+		{
+			n:      1,
+			dl:     nil,
+			d:      []float64{1.2, 2.3},
+			du:     nil,
+			panics: true,
+		},
+		{
+			n:      1,
+			dl:     []float64{},
+			d:      nil,
+			du:     []float64{},
+			panics: true,
+		},
+		{
+			n:      4,
+			dl:     []float64{1.2},
+			d:      nil,
+			du:     nil,
+			panics: true,
+		},
+		{
+			n:      4,
+			dl:     []float64{1.2, 2.3, 3.4},
+			d:      []float64{4.5, 5.6, 6.7, 7.8, 1.2},
+			du:     []float64{8.9, 9.0, 0.1},
+			panics: true,
+		},
+	} {
+		var a *Tridiag
+		panicked, msg := panics(func() {
+			a = NewTridiag(test.n, test.dl, test.d, test.du)
+		})
+		if panicked {
+			if !test.panics {
+				t.Errorf("Case %d: unexpected panic: %s", i, msg)
+			}
+			continue
+		}
+		if test.panics {
+			t.Errorf("Case %d: expected panic", i)
+			continue
+		}
+
+		r, c := a.Dims()
+		if r != test.n {
+			t.Errorf("Case %d: unexpected number of rows: got=%d want=%d", i, r, test.n)
+		}
+		if c != test.n {
+			t.Errorf("Case %d: unexpected number of columns: got=%d want=%d", i, c, test.n)
+		}
+
+		kl, ku := a.Bandwidth()
+		if kl != 1 || ku != 1 {
+			t.Errorf("Case %d: unexpected bandwidth: got=%d,%d want=1,1", i, kl, ku)
+		}
+
+		if !reflect.DeepEqual(a, test.want) {
+			t.Errorf("Case %d: unexpected value via reflect: got=%v, want=%v", i, a, test.want)
+		}
+		if !Equal(a, test.want) {
+			t.Errorf("Case %d: unexpected value via mat.Equal: got=%v, want=%v", i, a, test.want)
+		}
+		if !Equal(a, test.dense) {
+			t.Errorf("Case %d: unexpected value via mat.Equal(Tridiag,Dense):\ngot:\n% v\nwant:\n% v", i, Formatted(a), Formatted(test.dense))
+		}
+	}
+}
+
+func TestTridiagAtSet(t *testing.T) {
+	t.Parallel()
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		tri, ref := newTestTridiag(n)
+
+		name := fmt.Sprintf("Case n=%v", n)
+
+		// Check At explicitly with all valid indices.
+		for i := 0; i < n; i++ {
+			for j := 0; j < n; j++ {
+				if tri.At(i, j) != ref.At(i, j) {
+					t.Errorf("%v: unexpected value for At(%d,%d): got %v, want %v",
+						name, i, j, tri.At(i, j), ref.At(i, j))
+				}
+			}
+		}
+		// Check At via a call to Equal.
+		if !Equal(tri, ref) {
+			t.Errorf("%v: unexpected value:\ngot: % v\nwant:% v",
+				name, Formatted(tri, Prefix("     ")), Formatted(ref, Prefix("     ")))
+		}
+
+		// Check At out of bounds.
+		for _, i := range []int{-1, n, n + 1} {
+			for j := 0; j < n; j++ {
+				panicked, message := panics(func() { tri.At(i, j) })
+				if !panicked || message != ErrRowAccess.Error() {
+					t.Errorf("%v: expected panic for invalid row access at (%d,%d)", name, i, j)
+				}
+			}
+		}
+		for _, j := range []int{-1, n, n + 1} {
+			for i := 0; i < n; i++ {
+				panicked, message := panics(func() { tri.At(i, j) })
+				if !panicked || message != ErrColAccess.Error() {
+					t.Errorf("%v: expected panic for invalid column access at (%d,%d)", name, i, j)
+				}
+			}
+		}
+
+		// Check SetBand out of bounds.
+		for _, i := range []int{-1, n, n + 1} {
+			for j := 0; j < n; j++ {
+				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
+				if !panicked || message != ErrRowAccess.Error() {
+					t.Errorf("%v: expected panic for invalid row access at (%d,%d)", name, i, j)
+				}
+			}
+		}
+		for _, j := range []int{-1, n, n + 1} {
+			for i := 0; i < n; i++ {
+				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
+				if !panicked || message != ErrColAccess.Error() {
+					t.Errorf("%v: expected panic for invalid column access at (%d,%d)", name, i, j)
+				}
+			}
+		}
+		for i := 0; i < n; i++ {
+			for j := 0; j <= i-2; j++ {
+				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
+				if !panicked || message != ErrBandSet.Error() {
+					t.Errorf("%v: expected panic for invalid access at (%d,%d)", name, i, j)
+				}
+			}
+			for j := i + 2; j < n; j++ {
+				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
+				if !panicked || message != ErrBandSet.Error() {
+					t.Errorf("%v: expected panic for invalid access at (%d,%d)", name, i, j)
+				}
+			}
+		}
+
+		// Check SetBand within bandwidth.
+		for i := 0; i < n; i++ {
+			for j := max(0, i-1); j <= min(i+1, n-1); j++ {
+				want := float64(i*n + j + 100)
+				tri.SetBand(i, j, want)
+				if got := tri.At(i, j); got != want {
+					t.Errorf("%v: unexpected value at (%d,%d) after SetBand: got %v, want %v", name, i, j, got, want)
+				}
+			}
+		}
+	}
+}
+
+func newTestTridiag(n int) (*Tridiag, *Dense) {
+	var dl, d, du []float64
+	d = make([]float64, n)
+	if n > 1 {
+		dl = make([]float64, n-1)
+		du = make([]float64, n-1)
+	}
+	for i := range d {
+		d[i] = float64(i*n + i + 1)
+	}
+	for j := range dl {
+		i := j + 1
+		dl[j] = float64(i*n + j + 1)
+	}
+	for i := range du {
+		j := i + 1
+		du[i] = float64(i*n + j + 1)
+	}
+	dense := make([]float64, n*n)
+	for i := 0; i < n; i++ {
+		for j := max(0, i-1); j <= min(i+1, n-1); j++ {
+			dense[i*n+j] = float64(i*n + j + 1)
+		}
+	}
+	return NewTridiag(n, dl, d, du), NewDense(n, n, dense)
+}
+
+func TestTridiagReset(t *testing.T) {
+	t.Parallel()
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		a, _ := newTestTridiag(n)
+		if a.IsEmpty() {
+			t.Errorf("Case n=%d: matrix is empty", n)
+		}
+		a.Reset()
+		if !a.IsEmpty() {
+			t.Errorf("Case n=%d: matrix is not empty after Reset", n)
+		}
+	}
+}
+
+func TestTridiagDiagView(t *testing.T) {
+	t.Parallel()
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		a, _ := newTestTridiag(n)
+		testDiagView(t, n, a)
+	}
+}
+
+func TestTridiagZero(t *testing.T) {
+	t.Parallel()
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		a, _ := newTestTridiag(n)
+		a.Zero()
+		for i := 0; i < n; i++ {
+			for j := 0; j < n; j++ {
+				if a.At(i, j) != 0 {
+					t.Errorf("Case n=%d: unexpected non-zero at (%d,%d): got %f", n, i, j, a.At(i, j))
+				}
+			}
+		}
+	}
+}
+
+func TestTridiagTrace(t *testing.T) {
+	t.Parallel()
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		a, ref := newTestTridiag(n)
+		want := ref.Trace()
+		got := a.Trace()
+		if got != want {
+			t.Errorf("Case n=%d: unexpected trace: got %f, want %f", n, got, want)
+		}
+	}
+}
+
+func TestTridiagNorm(t *testing.T) {
+	t.Parallel()
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		for _, norm := range []float64{1, 2, math.Inf(1)} {
+			a, ref := newTestTridiag(n)
+			want := Norm(ref, norm)
+			got := a.Norm(norm)
+			if got != want {
+				t.Errorf("Case n=%d,norm=%.0f: unexpected norm: got %f, want %f", n, norm, got, want)
+			}
+		}
+	}
+}
+
+func TestTridiagSolveTo(t *testing.T) {
+	t.Parallel()
+
+	const tol = 1e-13
+
+	rnd := rand.New(rand.NewSource(1))
+	random := func(n int) []float64 {
+		d := make([]float64, n)
+		for i := range d {
+			d[i] = rnd.NormFloat64()
+		}
+		return d
+	}
+
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		a := NewTridiag(n, random(n-1), random(n), random(n-1))
+		var aDense Dense
+		aDense.CloneFrom(a)
+		for _, trans := range []bool{false, true} {
+			for _, nrhs := range []int{1, 2, 5} {
+				const (
+					denseB = iota
+					rawB
+					basicB
+				)
+				for _, bType := range []int{denseB, rawB, basicB} {
+					const (
+						emptyDst = iota
+						shapedDst
+						bIsDst
+					)
+					for _, dstType := range []int{emptyDst, shapedDst, bIsDst} {
+						if dstType == bIsDst && bType != denseB {
+							continue
+						}
+
+						var b Matrix
+						switch bType {
+						case denseB:
+							b = NewDense(n, nrhs, random(n*nrhs))
+						case rawB:
+							b = &rawMatrix{asBasicMatrix(NewDense(n, nrhs, random(n*nrhs)))}
+						case basicB:
+							b = asBasicMatrix(NewDense(n, nrhs, random(n*nrhs)))
+						default:
+							panic("bad bType")
+						}
+
+						var dst *Dense
+						switch dstType {
+						case emptyDst:
+							dst = new(Dense)
+						case shapedDst:
+							dst = NewDense(n, nrhs, random(n*nrhs))
+						case bIsDst:
+							dst = b.(*Dense)
+						default:
+							panic("bad dstType")
+						}
+
+						name := fmt.Sprintf("n=%d,nrhs=%d,trans=%t,dstType=%d,bType=%d", n, nrhs, trans, dstType, bType)
+
+						var want Dense
+						var err error
+						if !trans {
+							err = want.Solve(&aDense, b)
+						} else {
+							err = want.Solve(aDense.T(), b)
+						}
+						if err != nil {
+							t.Fatalf("%v: unexpected failure when computing reference solution: %v", name, err)
+						}
+
+						err = a.SolveTo(dst, trans, b)
+						if err != nil {
+							t.Fatalf("%v: unexpected failure from Tridiag.SolveTo: %v", name, err)
+						}
+
+						var diff Dense
+						diff.Sub(dst, &want)
+						if resid := Norm(&diff, 1); resid > tol*float64(n) {
+							t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol*float64(n))
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+func TestTridiagSolveVecTo(t *testing.T) {
+	t.Parallel()
+
+	const tol = 1e-13
+
+	rnd := rand.New(rand.NewSource(1))
+	random := func(n int) []float64 {
+		d := make([]float64, n)
+		for i := range d {
+			d[i] = rnd.NormFloat64()
+		}
+		return d
+	}
+
+	for _, n := range []int{1, 2, 3, 4, 7, 10} {
+		a := NewTridiag(n, random(n-1), random(n), random(n-1))
+		var aDense Dense
+		aDense.CloneFrom(a)
+		for _, trans := range []bool{false, true} {
+			const (
+				denseB = iota
+				rawB
+				basicB
+			)
+			for _, bType := range []int{denseB, rawB, basicB} {
+				const (
+					emptyDst = iota
+					shapedDst
+					bIsDst
+				)
+				for _, dstType := range []int{emptyDst, shapedDst, bIsDst} {
+					if dstType == bIsDst && bType != denseB {
+						continue
+					}
+
+					var b Vector
+					switch bType {
+					case denseB:
+						b = NewVecDense(n, random(n))
+					case rawB:
+						b = &rawVector{asBasicVector(NewVecDense(n, random(n)))}
+					case basicB:
+						b = asBasicVector(NewVecDense(n, random(n)))
+					default:
+						panic("bad bType")
+					}
+
+					var dst *VecDense
+					switch dstType {
+					case emptyDst:
+						dst = new(VecDense)
+					case shapedDst:
+						dst = NewVecDense(n, random(n))
+					case bIsDst:
+						dst = b.(*VecDense)
+					default:
+						panic("bad dstType")
+					}
+
+					name := fmt.Sprintf("n=%d,trans=%t,dstType=%d,bType=%d", n, trans, dstType, bType)
+
+					var want VecDense
+					var err error
+					if !trans {
+						err = want.SolveVec(&aDense, b)
+					} else {
+						err = want.SolveVec(aDense.T(), b)
+					}
+					if err != nil {
+						t.Fatalf("%v: unexpected failure when computing reference solution: %v", name, err)
+					}
+
+					err = a.SolveVecTo(dst, trans, b)
+					if err != nil {
+						t.Fatalf("%v: unexpected failure from Tridiag.SolveTo: %v", name, err)
+					}
+
+					var diff Dense
+					diff.Sub(dst, &want)
+					if resid := Norm(&diff, 1); resid > tol*float64(n) {
+						t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol*float64(n))
+					}
+				}
+			}
+		}
+	}
+}