diff --git a/cgo/lapack.go b/cgo/lapack.go
index ea5a30fd..37244390 100644
--- a/cgo/lapack.go
+++ b/cgo/lapack.go
@@ -20,6 +20,7 @@ const (
 	badDims         = "lapack: bad input dimensions"
 	badDirect       = "lapack: bad direct"
 	badE            = "lapack: e has insufficient length"
+	badEigComp      = "lapack: bad EigComp"
 	badIpiv         = "lapack: insufficient permutation length"
 	badLdA          = "lapack: index of a out of range"
 	badNorm         = "lapack: bad norm"
diff --git a/lapack.go b/lapack.go
index 0ae57efc..1d433129 100644
--- a/lapack.go
+++ b/lapack.go
@@ -109,3 +109,17 @@ const (
 	SVDOverwrite        = 'O' // Compute the singular vectors and store them in input matrix
 	SVDNone             = 'N' // Do not compute singular vectors
 )
+
+// EigComp specifies the type of eigenvalue decomposition.
+type EigComp byte
+
+const (
+	// EigValueOnly specifies to compute only the eigenvalues of the input matrix.
+	EigValueOnly EigComp = 'N'
+	// EigDecomp specifies to compute the eigenvalues and eigenvectors of the
+	// full symmetric matrix.
+	EigDecomp = 'V'
+	// EigBoth specifies to compute both the eigenvalues and eigenvectors of the
+	// input tridiagonal matrix.
+	EigBoth = 'I'
+)
diff --git a/native/dsteqr.go b/native/dsteqr.go
new file mode 100644
index 00000000..d408a7af
--- /dev/null
+++ b/native/dsteqr.go
@@ -0,0 +1,389 @@
+// Copyright ©2016 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 native
+
+import (
+	"math"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+	"github.com/gonum/lapack"
+)
+
+// Dsteqr computes the eigenvalues and optionally the eigenvectors of a symmetric
+// tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a
+// full or band symmetric matrix can also be found if Dsytrd, Dsptrd, or Dsbtrd
+// have been used to reduce this matrix to tridiagonal form.
+//
+// d, on entry, contains the diagonal elements of the tridiagonal matrix. On exit,
+// d contains the eigenvalues in ascending order. d must have length n and
+// Dsteqr will panic otherwise.
+//
+// e, on entry, contains the off-diagonal elements of the tridiagonal matrix on
+// entry, and is overwritten during the call to Dsteqr. e must have length n-1 and
+// Dsteqr will panic otherwise.
+//
+// z, on entry, contains the n×n orthogonal matrix used in the reduction to
+// tridiagonal form if compz == lapack.EigDecomp. On exit, if
+// compz == lapack.EigBoth, z contains the orthonormal eigenvectors of the
+// original symmetric matrix, and if compz == lapack.EigDecomp, z contains the
+// orthonormal eigenvectors of the symmetric tridiagonal matrix. z is not used
+// if compz == lapack.EigValueOnly.
+//
+// work must have length at least max(1, 2*n-2) if the eigenvectors are computed,
+// and Dsteqr will panic otherwise.
+func (impl Implementation) Dsteqr(compz lapack.EigComp, n int, d, e, z []float64, ldz int, work []float64) (ok bool) {
+	if len(d) < n {
+		panic(badD)
+	}
+	if len(e) < n-1 {
+		panic(badE)
+	}
+	if compz != lapack.EigValueOnly && compz != lapack.EigBoth && compz != lapack.EigDecomp {
+		panic(badEigComp)
+	}
+	if compz != lapack.EigValueOnly {
+		if len(work) < max(1, 2*n-2) {
+			panic(badWork)
+		}
+		checkMatrix(n, n, z, ldz)
+	}
+
+	var icompz int
+	if compz == lapack.EigDecomp {
+		icompz = 1
+	} else if compz == lapack.EigBoth {
+		icompz = 2
+	}
+
+	if n == 0 {
+		return true
+	}
+	if n == 1 {
+		if icompz == 2 {
+			z[0] = 1
+		}
+		return true
+	}
+
+	bi := blas64.Implementation()
+
+	eps := dlamchE
+	eps2 := eps * eps
+	safmin := dlamchS
+	safmax := 1 / safmin
+	ssfmax := math.Sqrt(safmax) / 3
+	ssfmin := math.Sqrt(safmin) / eps2
+
+	// Compute the eigenvalues and eigenvectors of the tridiagonal matrix.
+	if icompz == 2 {
+		impl.Dlaset(blas.All, n, n, 0, 1, z, ldz)
+	}
+	const maxit = 30
+	nmaxit := n * maxit
+
+	jtot := 0
+
+	// Determine where the matrix splits and choose QL or QR iteration for each
+	// block, according to whether top or bottom diagonal element is smaller.
+	// TODO(btracey): Reformat the gotos to use structural flow techniques.
+	// TODO(btracey): Move these variables closer to actual usage when
+	// gotos are removed.
+	l1 := 0
+	nm1 := n - 1
+	var l, lend, lendsv, lsv, m int
+	var anorm, c, g, p, r, s float64
+
+	type scaletype int
+	const (
+		none scaletype = iota
+		down
+		up
+	)
+	var iscale scaletype
+
+Ten:
+	if l1 > n-1 {
+		goto OneSixty
+	}
+	if l1 > 0 {
+		e[l1-1] = 0
+	}
+	if l1 <= nm1 {
+		for m = l1; m < nm1; m++ {
+			test := math.Abs(e[m])
+			if test == 0 {
+				goto Thirty
+			}
+			if test <= (math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1])))*eps {
+				e[m] = 0
+				goto Thirty
+			}
+		}
+	}
+	m = n - 1
+Thirty:
+	l = l1
+	lsv = l
+	lend = m
+	lendsv = lend
+	l1 = m + 1
+	if lend == l {
+		goto Ten
+	}
+
+	// Scale submatrix in rows and columns L to Lend
+	anorm = impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:])
+	switch {
+	case anorm == 0:
+		goto Ten
+	case anorm > ssfmax:
+		iscale = down
+		// TODO(btracey): Why is lda n?
+		impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], n)
+		impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], n)
+	case anorm < ssfmin:
+		iscale = up
+		// TODO(btracey): Why is lda n?
+		impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], n)
+		impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], n)
+	}
+
+	// Choose between QL and QR.
+	if math.Abs(d[lend]) < math.Abs(d[l]) {
+		lend = lsv
+		l = lendsv
+	}
+	if lend > l {
+		// QL Iteration. Look for small subdiagonal element.
+	Fourty:
+		if l != lend {
+			lendm1 := lend - 1
+			for m = l; m <= lendm1; m++ {
+				v := math.Abs(e[m])
+				if v*v <= (eps2*math.Abs(d[m]))*math.Abs(d[m+1])+safmin {
+					goto Sixty
+				}
+			}
+		}
+		m = lend
+	Sixty:
+		if m < lend {
+			e[m] = 0
+		}
+		p = d[l]
+		if m == l {
+			goto Eighty
+		}
+
+		// If remaining matrix is 2×2, use Dlae2 to compute its eigensystem.
+		if m == l+1 {
+			if icompz > 0 {
+				d[l], d[l+1], work[l], work[n-1+l] = impl.Dlaev2(d[l], e[l], d[l+1])
+				impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward,
+					n, 2, work[l:], work[n-1+l:], z[l:], ldz)
+			} else {
+				d[l], d[l+1] = impl.Dlae2(d[l], e[l], d[l+1])
+			}
+			e[l] = 0
+			l += 2
+			if l <= lend {
+				goto Fourty
+			}
+			goto OneFourty
+		}
+
+		if jtot == nmaxit {
+			goto OneFourty
+		}
+		jtot++
+
+		// Form shift
+		g = (d[l+1] - p) / (2 * e[l])
+		r = impl.Dlapy2(g, 1)
+		g = d[m] - p + e[l]/(g+math.Copysign(r, g))
+		s = 1
+		c = 1
+		p = 0
+
+		// Inner loop
+		for i := m - 1; i >= l; i-- {
+			f := s * e[i]
+			b := c * e[i]
+			c, s, r = impl.Dlartg(g, f)
+			if i != m-1 {
+				e[i+1] = r
+			}
+			g = d[i+1] - p
+			r = (d[i]-g)*s + 2*c*b
+			p = s * r
+			d[i+1] = g + p
+			g = c*r - b
+
+			// If eigenvectors are desired, then save rotations.
+			if icompz > 0 {
+				work[i] = c
+				work[n-1+i] = -s
+			}
+		}
+		// If eigenvectors are desired, then apply saved rotations.
+		if icompz > 0 {
+			mm := m - l + 1
+			impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward,
+				n, mm, work[l:], work[n-1+l:], z[l:], ldz)
+		}
+		d[l] -= p
+		e[l] = g
+		goto Fourty
+
+		// Eigenvalue found.
+	Eighty:
+		d[l] = p
+		l++
+		if l <= lend {
+			goto Fourty
+		}
+		goto OneFourty
+	} else {
+		// QR Iteration.
+		// Look for small superdiagonal element.
+	Ninety:
+		if l != lend {
+			lendp1 := lend + 1
+			for m = l; m >= lendp1; m-- {
+				v := math.Abs(e[m-1])
+				if v*v <= (eps2*math.Abs(d[m])*math.Abs(d[m-1]) + safmin) {
+					goto OneTen
+				}
+			}
+		}
+		m = lend
+	OneTen:
+		if m > lend {
+			e[m-1] = 0
+		}
+		p = d[l]
+		if m == l {
+			goto OneThirty
+		}
+
+		// If remaining matrix is 2×2, use Dlae2 to compute its eigenvalues.
+		if m == l-1 {
+			if icompz > 0 {
+				d[l-1], d[l], work[m], work[n-1+m] = impl.Dlaev2(d[l-1], e[l-1], d[l])
+				impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward,
+					n, 2, work[m:], work[n-1+m:], z[l-1:], ldz)
+			} else {
+				d[l-1], d[l] = impl.Dlae2(d[l-1], e[l-1], d[l])
+			}
+			e[l-1] = 0
+			l -= 2
+			if l >= lend {
+				goto Ninety
+			}
+			goto OneFourty
+		}
+		if jtot == nmaxit {
+			goto OneFourty
+		}
+		jtot++
+
+		// Form shift.
+		g = (d[l-1] - p) / (2 * e[l-1])
+		r = impl.Dlapy2(g, 1)
+		g = d[m] - p + (e[l-1])/(g+math.Copysign(r, g))
+
+		s = 1
+		c = 1
+		p = 0
+
+		// Inner loop.
+		for i := m; i < l; i++ {
+			f := s * e[i]
+			b := c * e[i]
+			c, s, r = impl.Dlartg(g, f)
+			if i != m {
+				e[i-1] = r
+			}
+			g = d[i] - p
+			r = (d[i+1]-g)*s + 2*c*b
+			p = s * r
+			d[i] = g + p
+			g = c*r - b
+
+			// If eigenvectors are desired, then save rotations.
+			if icompz > 0 {
+				work[i] = c
+				work[n-1+i] = s
+			}
+		}
+
+		// If eigenvectors are desired, then apply saved rotations.
+		if icompz > 0 {
+			mm := l - m + 1
+			impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward,
+				n, mm, work[m:], work[n-1+m:], z[m:], ldz)
+		}
+		d[l] -= p
+		e[l-1] = g
+		goto Ninety
+
+		// Eigenvalue found
+	OneThirty:
+		d[l] = p
+		l--
+		if l >= lend {
+			goto Ninety
+		}
+		goto OneFourty
+	}
+
+	// Undo scaling if necessary.
+OneFourty:
+	switch iscale {
+	case down:
+		impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], n)
+		impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv, 1, e[lsv:], n)
+	case up:
+		impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, e[lsv:], n)
+		impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv, 1, e[lsv:], n)
+	}
+
+	// Check for no convergence to an eigenvalue after a total of n*maxit iterations.
+	if jtot < nmaxit {
+		goto Ten
+	}
+	for i := 0; i < n-1; i++ {
+		if e[i] != 0 {
+			return false
+		}
+	}
+	return true
+
+	// Order eigenvalues and eigenvectors.
+OneSixty:
+	if icompz == 0 {
+		impl.Dlasrt(lapack.SortIncreasing, n, d)
+	} else {
+		// TODO(btracey): Consider replacing this sort with a call to sort.Sort.
+		for ii := 1; ii < n; ii++ {
+			i := ii - 1
+			k := i
+			p := d[i]
+			for j := ii; j < n; j++ {
+				if d[j] < p {
+					k = j
+					p = d[j]
+				}
+			}
+			if k != i {
+				d[k] = d[i]
+				d[i] = p
+				bi.Dswap(n, z[i:], ldz, z[k:], ldz)
+			}
+		}
+	}
+	return true
+}
diff --git a/native/general.go b/native/general.go
index eafd5465..a3bd2ded 100644
--- a/native/general.go
+++ b/native/general.go
@@ -26,6 +26,7 @@ const (
 	badDims         = "lapack: bad input dimensions"
 	badDirect       = "lapack: bad direct"
 	badE            = "lapack: e has insufficient length"
+	badEigComp      = "lapack: bad EigComp"
 	badIpiv         = "lapack: insufficient permutation length"
 	badLdA          = "lapack: index of a out of range"
 	badNorm         = "lapack: bad norm"
diff --git a/native/lapack_test.go b/native/lapack_test.go
index 1f535f65..24ab3279 100644
--- a/native/lapack_test.go
+++ b/native/lapack_test.go
@@ -228,6 +228,10 @@ func TestDrscl(t *testing.T) {
 	testlapack.DrsclTest(t, impl)
 }
 
+func TestDsteqr(t *testing.T) {
+	testlapack.DsteqrTest(t, impl)
+}
+
 func TestDsterf(t *testing.T) {
 	testlapack.DsterfTest(t, impl)
 }
diff --git a/testlapack/dsteqr.go b/testlapack/dsteqr.go
new file mode 100644
index 00000000..2f0bfdb6
--- /dev/null
+++ b/testlapack/dsteqr.go
@@ -0,0 +1,170 @@
+// Copyright ©2016 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 (
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+	"github.com/gonum/floats"
+	"github.com/gonum/lapack"
+)
+
+type Dsteqrer interface {
+	Dsteqr(compz lapack.EigComp, n int, d, e, z []float64, ldz int, work []float64) (ok bool)
+	Dorgtrer
+}
+
+func DsteqrTest(t *testing.T, impl Dsteqrer) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, compz := range []lapack.EigComp{lapack.EigDecomp, lapack.EigBoth} {
+		for _, test := range []struct {
+			n, lda int
+		}{
+			{1, 0},
+			{4, 0},
+			{8, 0},
+			{10, 0},
+
+			{2, 10},
+			{8, 10},
+			{10, 20},
+		} {
+			for cas := 0; cas < 100; cas++ {
+				n := test.n
+				lda := test.lda
+				if lda == 0 {
+					lda = n
+				}
+				d := make([]float64, n)
+				for i := range d {
+					d[i] = rnd.Float64()
+				}
+				e := make([]float64, n-1)
+				for i := range e {
+					e[i] = rnd.Float64()
+				}
+				a := make([]float64, n*lda)
+				for i := range a {
+					a[i] = rnd.Float64()
+				}
+				dCopy := make([]float64, len(d))
+				copy(dCopy, d)
+				eCopy := make([]float64, len(e))
+				copy(eCopy, e)
+				aCopy := make([]float64, len(a))
+				copy(aCopy, a)
+				if compz == lapack.EigDecomp {
+					// Compute triangular decomposition and orthonormal matrix.
+					uplo := blas.Upper
+					tau := make([]float64, n)
+					work := make([]float64, 1)
+					impl.Dsytrd(blas.Upper, n, a, lda, d, e, tau, work, -1)
+					work = make([]float64, int(work[0]))
+					impl.Dsytrd(uplo, n, a, lda, d, e, tau, work, len(work))
+					impl.Dorgtr(uplo, n, a, lda, tau, work, len(work))
+				} else {
+					for i := 0; i < n; i++ {
+						for j := 0; j < n; j++ {
+							a[i*lda+j] = 0
+							if i == j {
+								a[i*lda+j] = 1
+							}
+						}
+					}
+				}
+				work := make([]float64, 2*n)
+
+				aDecomp := make([]float64, len(a))
+				copy(aDecomp, a)
+				dDecomp := make([]float64, len(d))
+				copy(dDecomp, d)
+				eDecomp := make([]float64, len(e))
+				copy(eDecomp, e)
+				impl.Dsteqr(compz, n, d, e, a, lda, work)
+				dAns := make([]float64, len(d))
+				copy(dAns, d)
+
+				var truth blas64.General
+				if compz == lapack.EigDecomp {
+					truth = blas64.General{
+						Rows:   n,
+						Cols:   n,
+						Stride: n,
+						Data:   make([]float64, n*n),
+					}
+					for i := 0; i < n; i++ {
+						for j := i; j < n; j++ {
+							v := aCopy[i*lda+j]
+							truth.Data[i*truth.Stride+j] = v
+							truth.Data[j*truth.Stride+i] = v
+						}
+					}
+				} else {
+					truth = blas64.General{
+						Rows:   n,
+						Cols:   n,
+						Stride: n,
+						Data:   make([]float64, n*n),
+					}
+					for i := 0; i < n; i++ {
+						truth.Data[i*truth.Stride+i] = dCopy[i]
+						if i != n-1 {
+							truth.Data[(i+1)*truth.Stride+i] = eCopy[i]
+							truth.Data[i*truth.Stride+i+1] = eCopy[i]
+						}
+					}
+				}
+
+				V := blas64.General{
+					Rows:   n,
+					Cols:   n,
+					Stride: lda,
+					Data:   a,
+				}
+				if !eigenDecompCorrect(d, truth, V) {
+					t.Errorf("Eigen reconstruction mismatch. fromFull = %v, n = %v", compz == lapack.EigDecomp, n)
+				}
+
+				// Compare eigenvalues when not computing eigenvectors.
+				for i := range work {
+					work[i] = rnd.Float64()
+				}
+				impl.Dsteqr(lapack.EigValueOnly, n, dDecomp, eDecomp, aDecomp, lda, work)
+				if !floats.EqualApprox(d, dAns, 1e-8) {
+					t.Errorf("Eigenvalue mismatch when eigenvectors not computed")
+				}
+			}
+		}
+	}
+}
+
+// eigenDecompCorrect returns whether the eigen decomposition is correct.
+// It checks if
+//  A * v ≈ λ * v
+// where the eigenvalues λ are stored in values, and the eigenvectors are stored
+// in the columns of v.
+func eigenDecompCorrect(values []float64, A, V blas64.General) bool {
+	n := A.Rows
+	for i := 0; i < n; i++ {
+		lambda := values[i]
+		vector := make([]float64, n)
+		ans2 := make([]float64, n)
+		for j := range vector {
+			v := V.Data[j*V.Stride+i]
+			vector[j] = v
+			ans2[j] = lambda * v
+		}
+		v := blas64.Vector{Inc: 1, Data: vector}
+		ans1 := blas64.Vector{Inc: 1, Data: make([]float64, n)}
+		blas64.Gemv(blas.NoTrans, 1, A, v, 0, ans1)
+		if !floats.EqualApprox(ans1.Data, ans2, 1e-8) {
+			return false
+		}
+	}
+	return true
+}