diff --git a/lapack/gonum/dlag2.go b/lapack/gonum/dlag2.go
new file mode 100644
index 00000000..052dfb6d
--- /dev/null
+++ b/lapack/gonum/dlag2.go
@@ -0,0 +1,245 @@
+// 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"
+
+// Dlag2 computes the eigenvalues of a 2×2 generalized eigenvalue
+// problem
+//  A - w*B,
+// with scaling as necessary to avoid over-/underflow.
+// The scaling factor "s" results in a modified eigenvalue equation
+//  s*A - w*B
+// where  s  is a non-negative scaling factor chosen so that  w,  w*B,
+// and  s*A  do not overflow and, if possible, do not underflow, either.
+//
+// B is an upper triangular ldb×2 matrix
+// On entry, the 2×2 upper triangular matrix B.  It is
+// assumed that the one-norm of B is less than 1/dlamchS.  The
+// diagonals should be at least sqrt(dlamchS) times the largest
+// element of B (in absolute value); if a diagonal is smaller
+// than that, then  +/- sqrt(dlamchS) will be used instead of
+// that diagonal.
+//
+// It is assumed that A's 1-norm is less than 1/safmin.  Entries less than
+// sqrt(safmin)*norm(A) are subject to being treated as zero.
+//
+// Dlag2 is an internal routine. It is exported for testing purposes.
+//
+// scale1 is used to avoid over-/underflow in the
+// eigenvalue equation which defines the first eigenvalue.  If
+// the eigenvalues are complex, then the eigenvalues are
+// ( wr1  +/-  wi*i ) / scale1  (which may lie outside the
+// exponent range of the machine), scale1=scale2, and scale1
+// will always be positive.  If the eigenvalues are real, then
+// the first (real) eigenvalue is  wr1 / scale1 , but this may
+// overflow or underflow, and in fact, scale1 may be zero or
+// less than the underflow threshold if the exact eigenvalue
+// is sufficiently large.
+//
+// scale2 is used to avoid over-/underflow in the
+// eigenvalue equation which defines the second eigenvalue.  If
+// the eigenvalues are complex, then scale2=scale1.  If the
+// eigenvalues are real, then the second (real) eigenvalue is
+// wr2 / scale2 , but this may overflow or underflow, and in
+// fact, scale2 may be zero or less than the underflow
+// threshold if the exact eigenvalue is sufficiently large.
+//
+// If the eigenvalue is real, then wr1 is scale1 times the
+// eigenvalue closest to the (1,1) element of A*B^{-1}.  If the
+// eigenvalue is complex, then wr1=wr2 is scale1 times the real
+// part of the eigenvalues.
+//
+// If the eigenvalue is real, then wr2 is scale2 times the
+// other eigenvalue.  If the eigenvalue is complex, then
+// wr1=wr2 is scale1 times the real part of the eigenvalues.
+//
+// If the eigenvalue is real, then wi is zero.  If the
+// eigenvalue is complex, then wi is scale1 times the imaginary
+// part of the eigenvalues.  wi will always be non-negative.
+func (Implementation) Dlag2(a []float64, lda int, b []float64, ldb int) (scale1, scale2, wr1, wr2, wi float64) {
+
+	switch {
+	case lda < 2:
+		panic(badLdA)
+	case ldb < 2:
+		panic(badLdB)
+	case len(a) < lda+2:
+		panic(shortA)
+	case len(b) < ldb+2:
+		panic(shortB)
+	}
+
+	const (
+		safmin = dlamchS
+		fuzzy1 = 1 + 1e-5
+	)
+	rtmin := math.Sqrt(safmin)
+	rtmax := 1 / rtmin
+	safmax := 1 / safmin
+
+	// Scale a.
+	anorm := math.Max(math.Abs(a[0])+math.Abs(a[lda]),
+		math.Abs(a[1])+math.Abs(a[lda+1]))
+	anorm = math.Max(anorm, safmin)
+	ascale := 1 / anorm
+	a11 := ascale * a[0]
+	a21 := ascale * a[lda]
+	a12 := ascale * a[1]
+	a22 := ascale * a[lda+1]
+
+	// Perturb b if necessary to insure non-singularity.
+	b11 := b[0]
+	b12 := b[1]
+	b22 := b[ldb+1]
+	bmin := rtmin * math.Max(math.Max(math.Abs(b11), math.Abs(b12)),
+		math.Max(math.Abs(b22), rtmin))
+	if math.Abs(b11) < bmin {
+		b11 = math.Copysign(bmin, b11)
+	}
+	if math.Abs(b22) < bmin {
+		b22 = math.Copysign(bmin, b22)
+	}
+
+	// Scale B.
+	bnorm := math.Max(math.Max(math.Abs(b11), math.Abs(b12)+math.Abs(b22)), safmin)
+	bsize := math.Max(math.Abs(b11), math.Abs(b22))
+	bscale := 1 / bsize
+	b11 *= bscale
+	b12 *= bscale
+	b22 *= bscale
+
+	// Compute larger eigenvalue by method described by C. van Loan.
+	binv11 := 1 / b11
+	binv22 := 1 / b22
+	s1 := a11 * binv11
+	s2 := a22 * binv22
+	var as11, as12, as22, abi22, shift float64
+	var qq, ss, pp, discr, r float64
+	if math.Abs(s1) <= math.Abs(s2) {
+		as12 = a12 - s1*b12
+		as22 = a22 - s1*b22
+		ss = a21 * (binv11 * binv22)
+		abi22 = as22*binv22 - ss*b12
+		pp = 0.5 * abi22
+		shift = s1
+	} else {
+		as12 = a12 - s2*b12
+		as11 = a11 - s2*b11
+		ss = a21 * (binv11 * binv22)
+		abi22 = -ss * b12
+		pp = 0.5 * (as11*binv11 + abi22)
+		shift = s2
+	}
+	qq = ss * as12
+	pp2 := pp * pp
+	if math.Abs(pp*rtmin) >= 1 {
+		discr = rtmin*rtmin*pp2 + qq*safmin
+		r = math.Sqrt(math.Abs(discr)) * rtmax
+	} else {
+		if pp2+math.Abs(qq) <= safmin {
+			discr = rtmax*rtmax*pp2 + qq*safmax
+			r = math.Sqrt(math.Abs(discr)) * rtmin
+		} else {
+			discr = pp2 + qq
+			r = math.Sqrt(math.Abs(discr))
+		}
+	}
+
+	// Note: the test of R in the following `if` is to cover the case when
+	// discr is small and negative and is flushed to zero during
+	// the calculation of R.  On machines which have a consistent
+	// flush-to-zero threshold and handle numbers above that
+	// threshold correctly, it would not be necessary.
+	if discr >= 0 || r == 0 {
+		sum := pp + math.Copysign(r, pp)
+		diff := pp - math.Copysign(r, pp)
+		wbig := shift + sum
+
+		// Compute smaller eigenvalue.
+		wsmall := shift + diff
+		if 0.5*math.Abs(wbig) > math.Max(math.Abs(wsmall), safmin) {
+			wdet := (a11*a22 - a12*a21) * (binv11 * binv22)
+			wsmall = wdet / wbig
+		}
+		// Choose (real) eigenvalue closest to 2,2 element of A*B^{-1} for wr1.
+		if pp > abi22 {
+			wr1 = math.Min(wbig, wsmall)
+			wr2 = math.Max(wbig, wsmall)
+		} else {
+			wr1 = math.Max(wbig, wsmall)
+			wr2 = math.Min(wbig, wsmall)
+		}
+	} else {
+		// Complex eigenvalues.
+		wr1 = shift + pp
+		wr2 = wr1
+		wi = r
+	}
+
+	// Further scaling to avoid underflow and overflow in computing
+	// scale1 and overflow in computing w*B.
+	// This scale factor (wscale) is bounded from above using c1 and c2,
+	// and from below using c3 and c4.
+	// c1 implements the condition s*A must never overflow.
+	// c2 implements the condition w*B must never overflow.
+	// c3, with c2,
+	// implement the condition that s*A - w*B must never overflow.
+	// c4 implements the condition s should not underflow.
+	// c5 implements the condition max(s,|w|) should be at least 2.
+	c1 := bsize * (safmin * math.Max(1, ascale))
+	c2 := safmin * math.Max(1, bnorm)
+	c3 := bsize * safmin
+	c4 := 1.
+	c5 := 1.
+	if ascale <= 1 || bsize <= 1 {
+		c5 = math.Min(1, ascale*bsize)
+		if ascale <= 1 && bsize <= 1 {
+			c4 = math.Min(1, (ascale/safmin)*bsize)
+		}
+	}
+
+	// Scale first eigenvalue.
+	wabs := math.Abs(wr1) + math.Abs(wi)
+	wsize := math.Max(math.Max(safmin, c1), math.Max(fuzzy1*(wabs*c2+c3),
+		math.Min(c4, 0.5*math.Max(wabs, c5))))
+	maxABsize := math.Max(ascale, bsize)
+	minABsize := math.Min(ascale, bsize)
+	if wsize != 1 {
+		wscale := 1 / wsize
+		if wsize > 1 {
+			scale1 = (maxABsize * wscale) * minABsize
+		} else {
+			scale1 = (minABsize * wscale) * maxABsize
+		}
+		wr1 *= wscale
+		if wi != 0 {
+			wi *= wscale
+			wr2 = wr1
+			scale2 = scale1
+		}
+	} else {
+		scale1 = ascale * bsize
+		scale2 = scale1
+	}
+
+	// Scale second eigenvalue if real.
+	if wi == 0 {
+		wsize = math.Max(math.Max(safmin, c1), math.Max(fuzzy1*(math.Abs(wr2)*c2+c3),
+			math.Min(c4, 0.5*math.Max(math.Abs(wr2), c5))))
+		if wsize != 1 {
+			wscale := 1 / wsize
+			if wsize > 1 {
+				scale2 = (maxABsize * wscale) * minABsize
+			} else {
+				scale2 = (minABsize * wscale) * maxABsize
+			}
+			wr2 *= wscale
+		} else {
+			scale2 = ascale * bsize
+		}
+	}
+	return scale1, scale2, wr1, wr2, wi
+}
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index 5faa573b..8d30502e 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -193,6 +193,11 @@ func TestDlaexc(t *testing.T) {
 	testlapack.DlaexcTest(t, impl)
 }
 
+func TestDlag2(t *testing.T) {
+	t.Parallel()
+	testlapack.Dlag2Test(t, impl)
+}
+
 func TestDlags2(t *testing.T) {
 	t.Parallel()
 	testlapack.Dlags2Test(t, impl)
diff --git a/lapack/testlapack/dlag2.go b/lapack/testlapack/dlag2.go
new file mode 100644
index 00000000..0f402a8e
--- /dev/null
+++ b/lapack/testlapack/dlag2.go
@@ -0,0 +1,167 @@
+// 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/blas64"
+)
+
+type Dlag2er interface {
+	Dlag2(a []float64, lda int, b []float64, ldb int) (scale1, scale2, wr1, wr2, wi float64)
+}
+
+func Dlag2Test(t *testing.T, impl Dlag2er) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, lda := range []int{2, 5} {
+		for _, ldb := range []int{2, 5} {
+			for i := 0; i <= 80; i++ { // 80 iterations for variability.
+				dlag2Test(t, impl, rnd, lda, ldb)
+			}
+		}
+	}
+}
+
+func dlag2Test(t *testing.T, impl Dlag2er, rnd *rand.Rand, lda, ldb int) {
+	const tol = 1e-14
+	a := randomGeneral(2, 2, lda, rnd)
+	b := randomGeneral(2, 2, ldb, rnd)
+	a11 := a.Data[0]
+	a12 := a.Data[1]
+	a21 := a.Data[lda]
+	a22 := a.Data[lda+1]
+	b11 := b.Data[0]
+	b12 := b.Data[1]
+	b22 := b.Data[ldb+1]
+	name := fmt.Sprintf("lda=%d, ldb=%d", lda, ldb)
+	scale1, scale2, wr1, wr2, wi := impl.Dlag2(a.Data, a.Stride, b.Data, b.Stride)
+	if wi < 0 {
+		t.Errorf("%v: wi can not be negative wi=%v", name, wi)
+	}
+
+	dget53Test(t, a, b, scale1, wr1, wi)
+	dget53Test(t, a, b, scale2, wr2, wi)
+
+	// Complex eigenvalue solution.
+	if wi > 0 {
+		if wr1 != wr2 {
+			t.Errorf("%v: wr1 should equal wr2 for complex eigenvalues. got %v!=%v", name, wr1, wr2)
+		}
+		idet1 := cmplxdet2x2(complex(scale1*a11-wr1*b11, -wi*b11), complex(scale1*a12-wr1*b12, -wi*b12),
+			complex(scale1*a21, 0), complex(scale1*a22-wr1*b22, -wi*b22))
+		idet2 := cmplxdet2x2(complex(scale1*a11-wr1*b11, wi*b11), complex(scale1*a12-wr1*b12, wi*b12),
+			complex(scale1*a21, 0), complex(scale1*a22-wr1*b22, wi*b22))
+
+		if math.Abs(real(idet1))+math.Abs(imag(idet1)) > tol || math.Abs(real(idet2))+math.Abs(imag(idet2)) > tol {
+			t.Errorf("%v: (%.4v±%.4vi)/%.4v did not solve eigenvalue problem", name, wr1, wi, scale1)
+		}
+		return // Complex solution verified
+	}
+
+	// Real eigenvalue solution.
+	// |s*A - w*B| = 0   is solution for eigenvalue problem. Calculate distance from solution.
+	res1 := det2x2(scale1*a11-wr1*b11, scale1*a12-wr1*b12,
+		scale1*a21, scale1*a22-wr1*b22)
+	if math.Abs(res1) > tol {
+		t.Errorf("%v: got a residual from |%0.4v*A - w1*B| > tol: %v", name, scale1, res1)
+	}
+	res2 := det2x2(scale2*a11-wr2*b11, scale2*a12-wr2*b12,
+		scale2*a21, scale2*a22-wr2*b22)
+	if math.Abs(res2) > tol {
+		t.Errorf("%v: got a residual from |s*A - w2*B| > tol: %v", name, res1)
+	}
+}
+
+// dget53Test checks the generalized eigenvalues computed by dlag2.
+// The basic test for an eigenvalue is:
+//                           | det( s*A - w*B ) |
+//  result =  ---------------------------------------------------
+//            ulp max( s*norm(A), |w|*norm(B) )*norm( s*A - w*B )
+//
+// If s and w cant be scaled result will be 1/ulp.
+func dget53Test(t *testing.T, a, b blas64.General, scale, wr, wi float64) (result float64) {
+	lda := a.Stride
+	ldb := b.Stride
+	// Machine constants and norms.
+	safmin := dlamchS
+	ulp := dlamchE * dlamchB
+	absw := math.Abs(wr) + math.Abs(wi)
+	anorm := math.Max(math.Abs(a.Data[0])+math.Abs(a.Data[1]), math.Max(math.Abs(a.Data[lda])+math.Abs(a.Data[lda+1]), safmin))
+	bnorm := math.Max(math.Abs(b.Data[0]), math.Max(math.Abs(b.Data[ldb])+math.Abs(b.Data[ldb+1]), safmin))
+
+	// Check for possible overflow.
+	temp := (safmin*bnorm)*absw + (safmin*anorm)*scale
+
+	scales := scale
+	wrs := wr
+	wis := wi
+	if temp >= 1 {
+		t.Error("dget53: s*norm(A) + |w|*norm(B) > 1/safe_minimum")
+		temp = 1 / temp
+		scales *= temp
+		wrs *= temp
+		wis *= temp
+		absw = math.Abs(wrs) + math.Abs(wis)
+	}
+	s1 := math.Max(ulp*math.Max(scales*anorm, absw*bnorm), safmin*math.Max(scales, absw))
+
+	// Check for W and SCALE essentially zero.
+	if s1 < safmin {
+		t.Error("dget53: ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum")
+		if scales < safmin && absw < safmin {
+			t.Error("dget53: s and w could not be scaled so as to compute test")
+			return 1 / ulp
+		}
+		// Scale up to avoid underflow.
+		temp = 1 / math.Max(scales*anorm+absw*bnorm, safmin)
+		scales *= temp
+		wrs *= temp
+		wis *= temp
+		absw = math.Abs(wrs) + math.Abs(wis)
+		s1 = math.Max(ulp*math.Max(scales*anorm, absw*bnorm),
+			safmin*math.Max(scales, absw))
+		if s1 < safmin {
+			t.Error("dget53: s and w could not be scaled so as to compute test")
+			return 1 / ulp
+		}
+	}
+
+	// Compute C = s*A - w*B.
+	cr11 := scales*a.Data[0] - wrs*b.Data[0]
+	ci11 := -wis * b.Data[0]
+	cr21 := scales * a.Data[1]
+	cr12 := scales*a.Data[lda] - wrs*b.Data[ldb]
+	ci12 := -wis * b.Data[ldb]
+	cr22 := scales*a.Data[lda+1] - wrs*b.Data[ldb+1]
+	ci22 := -wis * b.Data[ldb+1]
+
+	// Compute the smallest singular value of s*A - w*B:
+	//                 |det( s*A - w*B )|
+	//     sigma_min = ------------------
+	//                 norm( s*A - w*B )
+	cnorm := math.Max(math.Abs(cr11)+math.Abs(ci11)+math.Abs(cr21),
+		math.Max(math.Abs(cr12)+math.Abs(ci12)+math.Abs(cr22)+math.Abs(ci22), safmin))
+	cscale := 1 / math.Sqrt(cnorm)
+	detr := (cscale*cr11)*(cscale*cr22) -
+		(cscale*ci11)*(cscale*ci22) -
+		(cscale*cr12)*(cscale*cr21)
+	deti := (cscale*cr11)*(cscale*ci22) +
+		(cscale*ci11)*(cscale*cr22) -
+		(cscale*ci12)*(cscale*cr21)
+	sigmin := math.Abs(detr) + math.Abs(deti)
+	result = sigmin / s1
+	return result
+}
+
+// cmplxdet2x2 returns the determinant of
+//  |a11 a12|
+//  |a21 a22|
+func cmplxdet2x2(a11, a12, a21, a22 complex128) complex128 {
+	return a11*a22 - a12*a21
+}