diff --git a/convex/lp/convert.go b/convex/lp/convert.go
new file mode 100644
index 0000000..ad897d0
--- /dev/null
+++ b/convex/lp/convert.go
@@ -0,0 +1,141 @@
+// 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 lp
+
+import (
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+// TODO(btracey): Have some sort of preprocessing step for helping to fix A to make it
+// full rank?
+// TODO(btracey): Reduce rows? Get rid of all zeros, places where only one variable
+// is there, etc. Could be implemented with a Reduce function.
+// TODO(btracey): Provide method of artificial variables for help when problem
+// is infeasible?
+// TODO(btracey): Add an lp.Solve that solves an LP in non-standard form.
+
+// Convert converts a General-form LP into a standard form LP.
+// The general form of an LP is:
+//  minimize c^T * x
+//  s.t      G * x <= h
+//           A * x = b
+// And the standard form is:
+//  minimize cNew^T * x
+//  s.t      aNew * x = bNew
+//           x >= 0
+// If there are no constraints of the given type, the inputs may be nil.
+func Convert(c []float64, g mat64.Matrix, h []float64, a mat64.Matrix, b []float64) (cNew []float64, aNew *mat64.Dense, bNew []float64) {
+	nVar := len(c)
+	nIneq := len(h)
+
+	// Check input sizes.
+	if g == nil {
+		if nIneq != 0 {
+			panic(badShape)
+		}
+	} else {
+		gr, gc := g.Dims()
+		if gr != nIneq {
+			panic(badShape)
+		}
+		if gc != nVar {
+			panic(badShape)
+		}
+	}
+
+	nEq := len(b)
+	if a == nil {
+		if nEq != 0 {
+			panic(badShape)
+		}
+	} else {
+		ar, ac := a.Dims()
+		if ar != nEq {
+			panic(badShape)
+		}
+		if ac != nVar {
+			panic(badShape)
+		}
+	}
+
+	// Convert the general form LP.
+	// Derivation:
+	// 0. Start with general form
+	//  min.	c^T * x
+	//  s.t.	G * x <= h
+	//  		A * x = b
+	// 1. Introduce slack variables for each constraint
+	//  min. 	c^T * x
+	//  s.t.	G * x + s = h
+	//			A * x = b
+	//      	s >= 0
+	// 2. Add non-negativity constraints for x by splitting x
+	// into positive and negative components.
+	//   x = xp - xn
+	//   xp >= 0, xn >= 0
+	// This makes the LP
+	//  min.	c^T * xp - c^T xn
+	//  s.t. 	G * xp - G * xn + s = h
+	//			A * xp  - A * xn = b
+	//			xp >= 0, xn >= 0, s >= 0
+	// 3. Write the above in standard form:
+	//  xt = [xp
+	//	 	  xn
+	//		  s ]
+	//  min.	[c^T, -c^T, 0] xt
+	//  s.t.	[G, -G, I] xt = h
+	//   		[A, -A, 0] xt = b
+	//			x >= 0
+
+	// In summary:
+	// Original LP:
+	//  min.	c^T * x
+	//  s.t.	G * x <= h
+	//  		A * x = b
+	// Standard Form:
+	//  xt = [xp; xn; s]
+	//  min.	[c^T, -c^T, 0] xt
+	//  s.t.	[G, -G, I] xt = h
+	//   		[A, -A, 0] xt = b
+	//			x >= 0
+
+	// New size of x is [xp, xn, s]
+	nNewVar := nVar + nVar + nIneq
+
+	// Construct cNew = [c; -c; 0]
+	cNew = make([]float64, nNewVar)
+	copy(cNew, c)
+	copy(cNew[nVar:], c)
+	floats.Scale(-1, cNew[nVar:2*nVar])
+
+	// New number of equality constraints is the number of total constraints.
+	nNewEq := nIneq + nEq
+
+	// Construct bNew = [h, b].
+	bNew = make([]float64, nNewEq)
+	copy(bNew, h)
+	copy(bNew[nIneq:], b)
+
+	// Construct aNew = [G, -G, I; A, -A, 0].
+	aNew = mat64.NewDense(nNewEq, nNewVar, nil)
+	if nIneq != 0 {
+		aView := (aNew.View(0, 0, nIneq, nVar)).(*mat64.Dense)
+		aView.Copy(g)
+		aView = (aNew.View(0, nVar, nIneq, nVar)).(*mat64.Dense)
+		aView.Scale(-1, g)
+		aView = (aNew.View(0, 2*nVar, nIneq, nIneq)).(*mat64.Dense)
+		for i := 0; i < nIneq; i++ {
+			aView.Set(i, i, 1)
+		}
+	}
+	if nEq != 0 {
+		aView := (aNew.View(nIneq, 0, nEq, nVar)).(*mat64.Dense)
+		aView.Copy(a)
+		aView = (aNew.View(nIneq, nVar, nEq, nVar)).(*mat64.Dense)
+		aView.Scale(-1, a)
+	}
+	return cNew, aNew, bNew
+}
diff --git a/convex/lp/simplex.go b/convex/lp/simplex.go
new file mode 100644
index 0000000..d8646f6
--- /dev/null
+++ b/convex/lp/simplex.go
@@ -0,0 +1,609 @@
+// 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 lp implements routines for solving linear programs.
+package lp
+
+import (
+	"errors"
+	"fmt"
+	"math"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+// TODO(btracey): Could have a solver structure with an abstract factorizer. With
+// this transformation the same high-level code could handle both Dense and Sparse.
+// TODO(btracey): Need to improve error handling. Only want to panic if condition number inf.
+// TODO(btracey): Performance enhancements. There are currently lots of linear
+// solves that can be improved by doing rank-one updates. For example, the swap
+// step is just a rank-one update.
+// TODO(btracey): Better handling on the linear solve errors. If the condition
+// number is not inf and the equation solved "well", should keep moving.
+
+var (
+	ErrBland      = errors.New("lp: bland: all replacements are negative or cause ill-conditioned ab")
+	ErrInfeasible = errors.New("lp: problem is infeasible")
+	ErrLinSolve   = errors.New("lp: linear solve failure")
+	ErrUnbounded  = errors.New("lp: problem is unbounded")
+	ErrSingular   = errors.New("lp: A is singular")
+	ErrZeroColumn = errors.New("lp: A has a column of all zeros")
+	ErrZeroRow    = errors.New("lp: A has a row of all zeros")
+)
+
+var (
+	badShape = "lp: size mismatch"
+)
+
+// TODO(btracey): Should these tolerances be part of a settings struct?
+
+const (
+	// initPosTol is the tolerance on the initial condition being feasible. Strictly,
+	// the x should be positive, but instead it must be greater than -initPosTol.
+	initPosTol = 1e-13
+	// blandNegTol is the tolerance on the value being greater than 0 in the bland test.
+	blandNegTol = 1e-14
+	// rRoundTol is the tolerance for rounding values to zero when testing if
+	// constraints are met.
+	rRoundTol = 1e-13
+	// dRoundTol is the tolerance for testing if values are zero for the problem
+	// being unbounded.
+	dRoundTol = 1e-13
+	// phaseIZeroTol tests if the Phase I problem returned a feasible solution.
+	phaseIZeroTol = 1e-12
+	// blandZeroTol is the tolerance on testing if the bland solution can move.
+	blandZeroTol = 1e-12
+)
+
+// Simplex solves a linear program in standard form using Danzig's Simplex
+// algorithm. The standard form of a linear program is:
+//  minimize	c^T x
+//  s.t. 		A*x = b
+//  			x >= 0 .
+// The input tol sets how close to the optimal solution is found (specifically,
+// when the maximal reduced cost is below tol). An error will be returned if the
+// problem is infeasible or unbounded. In rare cases, numeric errors can cause
+// the Simplex to fail. In this case, an error will be returned along with the
+// most recently found feasible solution.
+//
+// The Convert function can be used to transform a general LP into standard form.
+//
+// The input matrix A must have full rank and may not contain any columns with
+// all zeros. Furthermore, len(c) must equal the number of columns of A, and len(b)
+// must equal the number of rows of A. Simplex will panic if these conditions are
+// not met.
+//
+// initialBasic can be used to set the initial set of indices for a feasible
+// solution to the LP. If an initial feasible solution is not known, initialBasic
+// may be nil. If initialBasic is non-nil, len(initialBasic) must equal the number
+// of columns of A and must be an actual feasible solution to the LP, otherwise
+// Simplex will panic.
+//
+// A description of the Simplex algorithm can be found in Ch. 8 of
+//  Strang, Gilbert. "Linear Algebra and Applications." Academic, New York (1976).
+// For a detailed video introduction, see lectures 11-13 of UC Math 352
+//  https://www.youtube.com/watch?v=ESzYPFkY3og&index=11&list=PLh464gFUoJWOmBYla3zbZbc4nv2AXez6X.
+func Simplex(c []float64, A mat64.Matrix, b []float64, tol float64, initialBasic []int) (optF float64, optX []float64, err error) {
+	ans, x, _, err := simplex(initialBasic, c, A, b, tol)
+	return ans, x, err
+}
+
+func simplex(initialBasic []int, c []float64, A mat64.Matrix, b []float64, tol float64) (float64, []float64, []int, error) {
+	err := verifyInputs(initialBasic, c, A, b)
+	if err != nil {
+		if err == ErrUnbounded {
+			return math.Inf(-1), nil, nil, ErrUnbounded
+		}
+		return math.NaN(), nil, nil, err
+	}
+	m, n := A.Dims()
+
+	// There is at least one optimal solution to the LP which is at the intersection
+	// to a set of constraint boundaries. For a standard form LP with m variables
+	// and n equality constraints, at least m-n elements of x must equal zero
+	// at optimality. The Simplex algorithm solves the standard-form LP by starting
+	// at an initial constraint vertex and successively moving to adjacent constraint
+	// vertices. At every vertex, the set of non-zero x values is the "basic
+	// feasible solution". The list of non-zero x's are maintained in basicIdxs,
+	// the respective columns of A are in ab, and the actual non-zero values of
+	// x are in xb.
+	//
+	// The LP is equality constrained such that A * x = b. This can be expanded
+	// to
+	//  ab * xb + an * xn = b
+	// where ab are the columns of a in the basic set, and an are all of the
+	// other columns. Since each element of xn is zero by definition, this means
+	// that for all feasible solutions xb = ab^-1 * b.
+	//
+	// Before the simplex algorithm can start, an initial feasible solution must
+	// be found. If initialBasic is non-nil a feasible solution has been supplied.
+	// Otherwise the "Phase I" problem must be solved to find an initial feasible
+	// solution.
+
+	var basicIdxs []int // The indices of the non-zero x values.
+	var ab *mat64.Dense // The subset of columns of A listed in basicIdxs.
+	var xb []float64    // The non-zero elements of x. xb = ab^-1 b
+
+	if initialBasic != nil {
+		// InitialBasic supplied. Panic if incorrect length or infeasible.
+		if len(initialBasic) != m {
+			panic("lp: incorrect number of initial vectors")
+		}
+		ab = extractColumns(A, initialBasic)
+		xb, err = initializeFromBasic(ab, b)
+		if err != nil {
+			panic(err)
+		}
+		basicIdxs = make([]int, len(initialBasic))
+		copy(basicIdxs, initialBasic)
+	} else {
+		// No inital basis supplied. Solve the PhaseI problem.
+		basicIdxs, ab, xb, err = findInitialBasic(A, b)
+		if err != nil {
+			return math.NaN(), nil, nil, err
+		}
+	}
+
+	// basicIdxs contains the indexes for an initial feasible solution,
+	// ab contains the extracted columns of A, and xb contains the feasible
+	// solution. All x not in the basic set are 0 by construction.
+
+	// nonBasicIdx is the set of nonbasic variables.
+	nonBasicIdx := make([]int, 0, n-m)
+	inBasic := make(map[int]struct{})
+	for _, v := range basicIdxs {
+		inBasic[v] = struct{}{}
+	}
+	for i := 0; i < n; i++ {
+		_, ok := inBasic[i]
+		if !ok {
+			nonBasicIdx = append(nonBasicIdx, i)
+		}
+	}
+
+	// cb is the subset of c for the basic variables. an and cn
+	// are the equivalents to ab and cb but for the nonbasic variables.
+	cb := make([]float64, len(basicIdxs))
+	for i, idx := range basicIdxs {
+		cb[i] = c[idx]
+	}
+	cn := make([]float64, len(nonBasicIdx))
+	for i, idx := range nonBasicIdx {
+		cn[i] = c[idx]
+	}
+	an := extractColumns(A, nonBasicIdx)
+
+	bVec := mat64.NewVector(len(b), b)
+	cbVec := mat64.NewVector(len(cb), cb)
+
+	// Temporary data needed each iteration. (Described later)
+	r := make([]float64, n-m)
+	move := make([]float64, m)
+
+	// Solve the linear program starting from the initial feasible set. This is
+	// the "Phase 2" problem.
+	//
+	// Algorithm:
+	// 1) Compute the "reduced costs" for the non-basic variables. The reduced
+	// costs are the lagrange multipliers of the constraints.
+	// 	 r = cn - an^T * ab^-T * cb
+	// 2) If all of the reduced costs are positive, no improvement is possible,
+	// and the solution is optimal (xn can only increase because of
+	// non-negativity constraints). Otherwise, the solution can be improved and
+	// one element will be exchanged in the basic set.
+	// 3) Choose the x_n with the most negative value of r. Call this value xe.
+	// This variable will be swapped into the basic set.
+	// 4) Increase xe until the next constraint boundary is met. This will happen
+	// when the first element in xb becomes 0. The distance xe can increase before
+	// a given element in xb becomes negative can be found from
+	//	xb = Ab^-1 b - Ab^-1 An xn
+	//     = Ab^-1 b - Ab^-1 Ae xe
+	//     = bhat + d x_e
+	//  xe = bhat_i / - d_i
+	// where Ae is the column of A corresponding to xe.
+	// The constraining basic index is the first index for which this is true,
+	// so remove the element which is min_i (bhat_i / -d_i), assuming d_i is negative.
+	// If no d_i is less than 0, then the problem is unbounded.
+	// 5) If the new xe is 0 (that is, bhat_i == 0), then this location is at
+	// the intersection of several constraints. Use the Bland rule instead
+	// of the rule in step 4 to avoid cycling.
+	for {
+		// Compute reduced costs -- r = cn - an^T ab^-T cb
+		var tmp mat64.Vector
+		err = tmp.SolveVec(ab.T(), cbVec)
+		if err != nil {
+			break
+		}
+		data := make([]float64, n-m)
+		tmp2 := mat64.NewVector(n-m, data)
+		tmp2.MulVec(an.T(), &tmp)
+		floats.SubTo(r, cn, data)
+
+		// Replace the most negative element in the simplex. If there are no
+		// negative entries then the optimal solution has been found.
+		minIdx := floats.MinIdx(r)
+		if r[minIdx] >= -tol {
+			break
+		}
+
+		for i, v := range r {
+			if math.Abs(v) < rRoundTol {
+				r[i] = 0
+			}
+		}
+
+		// Compute the moving distance.
+		err = computeMove(move, minIdx, A, ab, xb, nonBasicIdx)
+		if err != nil {
+			if err == ErrUnbounded {
+				return math.Inf(-1), nil, nil, ErrUnbounded
+			}
+			break
+		}
+
+		// Replace the basic index along the tightest constraint.
+		replace := floats.MinIdx(move)
+		if move[replace] <= 0 {
+			replace, minIdx, err = replaceBland(A, ab, xb, basicIdxs, nonBasicIdx, r, move)
+			if err != nil {
+				if err == ErrUnbounded {
+					return math.Inf(-1), nil, nil, ErrUnbounded
+				}
+				break
+			}
+		}
+
+		// Replace the constrained basicIdx with the newIdx.
+		basicIdxs[replace], nonBasicIdx[minIdx] = nonBasicIdx[minIdx], basicIdxs[replace]
+		cb[replace], cn[minIdx] = cn[minIdx], cb[replace]
+		tmpCol1 := mat64.Col(nil, replace, ab)
+		tmpCol2 := mat64.Col(nil, minIdx, an)
+		ab.SetCol(replace, tmpCol2)
+		an.SetCol(minIdx, tmpCol1)
+
+		// Compute the new xb.
+		xbVec := mat64.NewVector(len(xb), xb)
+		err = xbVec.SolveVec(ab, bVec)
+		if err != nil {
+			break
+		}
+	}
+	// Found the optimum successfully or died trying. The basic variables get
+	// their values, and the non-basic variables are all zero.
+	opt := floats.Dot(cb, xb)
+	xopt := make([]float64, n)
+	for i, v := range basicIdxs {
+		xopt[v] = xb[i]
+	}
+	return opt, xopt, basicIdxs, err
+}
+
+// computeMove computes how far can be moved replacing each index. The results
+// are stored into move.
+func computeMove(move []float64, minIdx int, A mat64.Matrix, ab *mat64.Dense, xb []float64, nonBasicIdx []int) error {
+	// Find ae.
+	col := mat64.Col(nil, nonBasicIdx[minIdx], A)
+	aCol := mat64.NewVector(len(col), col)
+
+	// d = - Ab^-1 Ae
+	nb, _ := ab.Dims()
+	d := make([]float64, nb)
+	dVec := mat64.NewVector(nb, d)
+	err := dVec.SolveVec(ab, aCol)
+	if err != nil {
+		return ErrLinSolve
+	}
+	floats.Scale(-1, d)
+
+	for i, v := range d {
+		if math.Abs(v) < dRoundTol {
+			d[i] = 0
+		}
+	}
+
+	// If no di < 0, then problem is unbounded.
+	if floats.Min(d) >= 0 {
+		return ErrUnbounded
+	}
+
+	// move = bhat_i / - d_i, assuming d is negative.
+	bHat := xb // ab^-1 b
+	for i, v := range d {
+		if v >= 0 {
+			move[i] = math.Inf(1)
+		} else {
+			move[i] = bHat[i] / math.Abs(v)
+		}
+	}
+	return nil
+}
+
+// replaceBland uses the Bland rule to find the indices to swap if the minimum
+// move is 0. The indices to be swapped are replace and minIdx (following the
+// nomenclature in the main routine).
+func replaceBland(A mat64.Matrix, ab *mat64.Dense, xb []float64, basicIdxs, nonBasicIdx []int, r, move []float64) (replace, minIdx int, err error) {
+	// Use the traditional bland rule, except don't replace a constraint which
+	// causes the new ab to be singular.
+	for i, v := range r {
+		if v > -blandNegTol {
+			continue
+		}
+		minIdx = i
+		err = computeMove(move, minIdx, A, ab, xb, nonBasicIdx)
+		if err != nil {
+			// Either unbounded or something went wrong.
+			return -1, -1, err
+		}
+		replace = floats.MinIdx(move)
+		if math.Abs(move[replace]) > blandZeroTol {
+			// Large enough that it shouldn't be a problem
+			return replace, minIdx, nil
+		}
+		// Find a zero index where replacement is non-singular.
+		biCopy := make([]int, len(basicIdxs))
+		for replace, v := range move {
+			if v > blandZeroTol {
+				continue
+			}
+			copy(biCopy, basicIdxs)
+			biCopy[replace] = nonBasicIdx[minIdx]
+			abTmp := extractColumns(A, biCopy)
+			// If the condition number is reasonable, use this index.
+			if mat64.Cond(abTmp, 1) < 1e16 {
+				return replace, minIdx, nil
+			}
+		}
+	}
+	return -1, -1, ErrBland
+}
+
+func verifyInputs(initialBasic []int, c []float64, A mat64.Matrix, b []float64) error {
+	m, n := A.Dims()
+	if len(c) != n {
+		panic("lp: c vector incorrect length")
+	}
+	if len(b) != m {
+		panic("lp: b vector incorrect length")
+	}
+	if len(c) != n {
+		panic("lp: c vector incorrect length")
+	}
+	if len(initialBasic) != 0 && len(initialBasic) != m {
+		panic("lp: initialBasic incorrect length")
+	}
+
+	// Do some sanity checks so that ab does not become singular during the
+	// simplex solution. If the ZeroRow checks are removed then the code for
+	// finding a set of linearly indepent columns must be improved.
+
+	// Check that if a row of A only has zero elements that corresponding
+	// element in b is zero, otherwise the problem is infeasible.
+	// Otherwise return ErrZeroRow.
+	for i := 0; i < m; i++ {
+		isZero := true
+		for j := 0; j < n; j++ {
+			if A.At(i, j) != 0 {
+				isZero = false
+				break
+			}
+		}
+		if isZero && b[i] != 0 {
+			// Infeasible
+			return ErrInfeasible
+		} else if isZero {
+			return ErrZeroRow
+		}
+	}
+	// Check that if a column only has zero elements that the respective C vector
+	// is positive (otherwise unbounded). Otherwise return ErrZeroColumn.
+	for j := 0; j < n; j++ {
+		isZero := true
+		for i := 0; i < m; i++ {
+			if A.At(i, j) != 0 {
+				isZero = false
+				break
+			}
+		}
+		if isZero && c[j] < 0 {
+			return ErrUnbounded
+		} else if isZero {
+			return ErrZeroColumn
+		}
+	}
+	return nil
+}
+
+// initializeFromBasic initializes the basic feasible solution given a set of
+// basic indices. It extracts the columns of A specified by basicIdxs and finds
+// the x values at that location. If the columns of A are not linearly independent
+// or if the initial set is not feasible, valid is false.
+func initializeFromBasic(ab *mat64.Dense, b []float64) (xb []float64, err error) {
+	m, _ := ab.Dims()
+	xb = make([]float64, m)
+	xbMat := mat64.NewVector(m, xb)
+
+	err = xbMat.SolveVec(ab, mat64.NewVector(m, b))
+	if err != nil {
+		return nil, errors.New("lp: subcolumns of A for supplied initial basic singular")
+	}
+	// The solve ensures that the equality constraints are met (ab * xb = b).
+	// Thus, the solution is feasible if and only if all of the x's are positive.
+	allPos := true
+	for _, v := range xb {
+		if v < -initPosTol {
+			allPos = false
+			break
+		}
+	}
+	if !allPos {
+		return xb, errors.New("lp: supplied subcolumns not a feasible solution")
+	}
+	return xb, nil
+}
+
+// extractColumns creates a new matrix out of the columns of A specified by cols.
+// TODO(btracey): Allow this to take a receiver.
+func extractColumns(A mat64.Matrix, cols []int) *mat64.Dense {
+	r, _ := A.Dims()
+	sub := mat64.NewDense(r, len(cols), nil)
+	col := make([]float64, r)
+	for j, idx := range cols {
+		mat64.Col(col, idx, A)
+		sub.SetCol(j, col)
+	}
+	return sub
+}
+
+// findInitialBasic finds an initial basic solution, and returns the basic
+// indices, ab, and xb.
+func findInitialBasic(A mat64.Matrix, b []float64) ([]int, *mat64.Dense, []float64, error) {
+	m, n := A.Dims()
+	basicIdxs := findLinearlyIndependent(A)
+	if len(basicIdxs) != m {
+		return nil, nil, nil, ErrSingular
+	}
+
+	// It may be that this linearly independent basis is also a feasible set. If
+	// so, the Phase I problem can be avoided.
+	ab := extractColumns(A, basicIdxs)
+	xb, err := initializeFromBasic(ab, b)
+	if err == nil {
+		return basicIdxs, ab, xb, nil
+	}
+
+	// This set was not feasible. Instead the "Phase I" problem must be solved
+	// to find an initial feasible set of basis.
+	//
+	// Method: Construct an LP whose optimal solution is a feasible solution
+	// to the original LP.
+	// 1) Introduce an artificial variable x_{n+1}.
+	// 2) Let x_j be the most negative element of x_b (largest constraint violation).
+	// 3) Add the artificial variable to A with:
+	//      a_{n+1} = b - \sum_{i in basicIdxs} a_i + a_j
+	//    swap j with n+1 in the basicIdxs.
+	// 4) Define a new LP:
+	//   minimize  x_{n+1}
+	//   subject to [A A_{n+1}][x_1 ... x_{n+1}] = b
+	//          x, x_{n+1} >= 0
+	// 5) Solve this LP. If x_{n+1} != 0, then the problem is infeasible, otherwise
+	// the found basis can be used as an initial basis for phase II.
+	//
+	// The extra column in Step 3 is defined such that the vector of 1s is an
+	// initial feasible solution.
+
+	// Find the largest constraint violator.
+	// Compute a_{n+1} = b - \sum{i in basicIdxs}a_i + a_j. j is in basicIDx, so
+	// instead just subtract the basicIdx columns that are not minIDx.
+	minIdx := floats.MinIdx(xb)
+	aX1 := make([]float64, m)
+	copy(aX1, b)
+	col := make([]float64, m)
+	for i, v := range basicIdxs {
+		if i == minIdx {
+			continue
+		}
+		mat64.Col(col, v, A)
+		floats.Sub(aX1, col)
+	}
+
+	// Construct the new LP.
+	// aNew = [A, a_{n+1}]
+	// bNew = b
+	// cNew = 1 for x_{n+1}
+	aNew := mat64.NewDense(m, n+1, nil)
+	aNew.Copy(A)
+	aNew.SetCol(n, aX1)
+	basicIdxs[minIdx] = n // swap minIdx with n in the basic set.
+	c := make([]float64, n+1)
+	c[n] = 1
+
+	// Solve the Phase 2 linear program.
+	_, xOpt, newBasic, err := simplex(basicIdxs, c, aNew, b, 1e-10)
+	if err != nil {
+		return nil, nil, nil, errors.New(fmt.Sprintf("lp: error finding feasible basis: %s", err))
+	}
+
+	// If n+1 is part of the solution basis then the problem is infeasible. If
+	// not, then the problem is feasible and newBasic is an initial feasible
+	// solution.
+	if math.Abs(xOpt[n]) > phaseIZeroTol {
+		return nil, nil, nil, ErrInfeasible
+	}
+
+	// The value is zero. First, see if it's not in the basis (feasible solution).
+	basicIdx := -1
+	basicMap := make(map[int]struct{})
+	for i, v := range newBasic {
+		if v == n {
+			basicIdx = i
+		}
+		basicMap[v] = struct{}{}
+		xb[i] = xOpt[v]
+	}
+	if basicIdx == -1 {
+		// Not in the basis.
+		ab = extractColumns(A, newBasic)
+		return newBasic, ab, xb, nil
+	}
+
+	// The value is zero, but it's in the basis. See if swapping in another column
+	// finds a feasible solution.
+	for i := range xOpt {
+		if _, inBasic := basicMap[i]; inBasic {
+			continue
+		}
+		newBasic[basicIdx] = i
+		ab := extractColumns(A, newBasic)
+		xb, err := initializeFromBasic(ab, b)
+		if err == nil {
+			return newBasic, ab, xb, nil
+		}
+	}
+	return nil, nil, nil, ErrInfeasible
+}
+
+// findLinearlyIndependnt finds a set of linearly independent columns of A, and
+// returns the column indexes of the linearly independent columns.
+func findLinearlyIndependent(A mat64.Matrix) []int {
+	m, n := A.Dims()
+	idxs := make([]int, 0, m)
+	columns := mat64.NewDense(m, m, nil)
+	newCol := make([]float64, m)
+	// Walk in reverse order because slack variables are typically the last columns
+	// of A.
+	for i := n - 1; i >= 0; i-- {
+		if len(idxs) == m {
+			break
+		}
+		mat64.Col(newCol, i, A)
+		if len(idxs) == 0 {
+			// A column is linearly independent from the null set.
+			// This is what needs to be changed if zero columns are allowed, as
+			// a column of all zeros is not linearly independent from itself.
+			columns.SetCol(len(idxs), newCol)
+			idxs = append(idxs, i)
+			continue
+		}
+		if linearlyDependent(mat64.NewVector(m, newCol), columns.View(0, 0, m, len(idxs))) {
+			continue
+		}
+		columns.SetCol(len(idxs), newCol)
+		idxs = append(idxs, i)
+	}
+	return idxs
+}
+
+// linearlyDependent returns whether the vector is linearly dependent
+// with the columns of A. It assumes that A is a full-rank matrix.
+func linearlyDependent(vec *mat64.Vector, A mat64.Matrix) bool {
+	// Add vec to the columns of A, and see if the condition number is reasonable.
+	m, n := A.Dims()
+	aNew := mat64.NewDense(m, n+1, nil)
+	aNew.Copy(A)
+	col := mat64.Col(nil, 0, vec)
+	aNew.SetCol(n, col)
+	cond := mat64.Cond(aNew, 1)
+	return cond > 1e12
+}
diff --git a/convex/lp/simplex_test.go b/convex/lp/simplex_test.go
new file mode 100644
index 0000000..e5e1b66
--- /dev/null
+++ b/convex/lp/simplex_test.go
@@ -0,0 +1,272 @@
+// 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 lp
+
+import (
+	"fmt"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+const convergenceTol = 1e-10
+
+func TestSimplex(t *testing.T) {
+	// First test specific inputs. These were collected from failures
+	// during randomized testing.
+	// TODO(btracey): Test specific problems with known solutions.
+	for _, test := range []struct {
+		A            mat64.Matrix
+		b            []float64
+		c            []float64
+		tol          float64
+		initialBasic []int
+	}{
+		{
+			// Basic feasible LP
+			A: mat64.NewDense(2, 4, []float64{
+				-1, 2, 1, 0,
+				3, 1, 0, 1,
+			}),
+			b: []float64{4, 9},
+			c: []float64{-1, -2, 0, 0},
+			//initialBasic: nil,
+			tol: 0,
+		},
+		{
+			// Zero row that caused linear solver failure
+			A: mat64.NewDense(3, 5, []float64{0.09917822373225804, 0, 0, -0.2588175087223661, -0.5935518220870567, 1.301111422556007, 0.12220247487326946, 0, 0, -1.9194869979254463, 0, 0, 0, 0, -0.8588221231396473}),
+			b: []float64{0, 0, 0},
+			c: []float64{0, 0.598992624019304, 0, 0, 0},
+		},
+		{
+			// Case that caused linear solver failure
+			A: mat64.NewDense(13, 26, []float64{-0.7001209024399848, -0.7027502615621812, -0, 0.7354444798695736, -0, 0.476457578966189, 0.7001209024399848, 0.7027502615621812, 0, -0.7354444798695736, 0, -0.476457578966189, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.8446087238391438, -0, -0, 0.7705609478497938, -0, -0, -2.7311218710244463, 0, 0, -0.7705609478497938, 0, 0, 2.7311218710244463, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -0, 0.8332519091897401, 0.7762132098737671, -0, -0, -0.052470638647269585, 0, -0.8332519091897401, -0.7762132098737671, 0, 0, 0.052470638647269585, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.9898577208292023, 0.31653724289408824, -0, -0, 0.17797227766447388, 1.2702427184954932, -0.7998764021535656, -0.31653724289408824, 0, 0, -0.17797227766447388, -1.2702427184954932, 0.7998764021535656, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -0, -0, -0, -0, 0.4206278126213235, -0.7253374879437113, 0, 0, 0, 0, -0.4206278126213235, 0.7253374879437113, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, -0, -0.7567988418466963, 0.3304567624749696, 0.8385927625193501, -0.0021606686026376387, -0, 0, 0.7567988418466963, -0.3304567624749696, -0.8385927625193501, 0.0021606686026376387, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -2.230107839590404, -0.9897104202085316, -0, 0.24703471683023603, -0, -2.382860345431941, 0.6206871648345162, 0.9897104202085316, 0, -0.24703471683023603, 0, 2.382860345431941, -0.6206871648345162, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 1.4350469221322282, -0.9730343818431852, -0, 2.326429855201535, -0, -0.14347849887004038, -1.4350469221322282, 0.9730343818431852, 0, -2.326429855201535, 0, 0.14347849887004038, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -0.7943912888763849, -0.13735037357335078, -0.5101161104860161, -0, -0, -1.4790634590370297, 0.050911195996747316, 0.13735037357335078, 0.5101161104860161, 0, 0, 1.4790634590370297, -0.050911195996747316, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, -0, -0.2515400440591492, 0.2058339272568599, -0, -0, -1.314023802253438, 0, 0.2515400440591492, -0.2058339272568599, 0, 0, 1.314023802253438, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0.08279503413614919, -0.16669071891829756, -0, -0.6208413721884664, -0, -0.6348258970402827, -0.08279503413614919, 0.16669071891829756, 0, 0.6208413721884664, 0, 0.6348258970402827, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, -0, -0, 0.49634739711260845, -0, -0, -0, 0, 0, -0.49634739711260845, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1, -0.2797437186631715, -0.8356683570259136, 1.8970426594969672, -0.4095711945594497, 0.45831284820623924, -0.6109615338552246, 0.2797437186631715, 0.8356683570259136, -1.8970426594969672, 0.4095711945594497, -0.45831284820623924, 0.6109615338552246, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1}),
+			b: []float64{-0.8446087238391436, 0, 1.9898577208292023, 0, 0, -2.230107839590404, 0, 0.20560871112361512, 0, 0, 0, 0, 0},
+			c: []float64{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},
+		},
+		{
+			// Phase 1 of the above that panicked
+			A: mat64.NewDense(26, 52, []float64{0.7001209024399848, -0, -0, -0.31653724289408824, -0, -0, 0.9897104202085316, -1.4350469221322282, 0.13735037357335078, -0, -0.08279503413614919, -0, 0.2797437186631715, -0.7001209024399848, 0, 0, 0.31653724289408824, 0, 0, -0.9897104202085316, 1.4350469221322282, -0.13735037357335078, 0, 0.08279503413614919, 0, -0.2797437186631715, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7027502615621812, -0, -0.8332519091897401, -0, -0, 0.7567988418466963, -0, 0.9730343818431852, 0.5101161104860161, 0.2515400440591492, 0.16669071891829756, -0, 0.8356683570259136, -0.7027502615621812, 0, 0.8332519091897401, 0, 0, -0.7567988418466963, 0, -0.9730343818431852, -0.5101161104860161, -0.2515400440591492, -0.16669071891829756, 0, -0.8356683570259136, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0.7705609478497938, -0.7762132098737671, -0, -0, -0.3304567624749696, -0.24703471683023603, -0, -0, -0.2058339272568599, -0, -0.49634739711260845, -1.8970426594969672, 0, 0.7705609478497938, 0.7762132098737671, 0, 0, 0.3304567624749696, 0.24703471683023603, 0, 0, 0.2058339272568599, 0, 0.49634739711260845, 1.8970426594969672, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7354444798695736, -0, -0, -0.17797227766447388, -0, -0.8385927625193501, -0, -2.326429855201535, -0, -0, 0.6208413721884664, -0, 0.4095711945594497, 0.7354444798695736, 0, 0, 0.17797227766447388, 0, 0.8385927625193501, 0, 2.326429855201535, 0, 0, -0.6208413721884664, 0, -0.4095711945594497, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -1.2702427184954932, -0.4206278126213235, 0.0021606686026376387, 2.382860345431941, -0, 1.4790634590370297, -0, -0, -0, -0.45831284820623924, 0, 0, 0, 1.2702427184954932, 0.4206278126213235, -0.0021606686026376387, -2.382860345431941, 0, -1.4790634590370297, 0, 0, 0, 0.45831284820623924, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.476457578966189, 2.7311218710244463, 0.052470638647269585, 0.7998764021535656, 0.7253374879437113, -0, -0.6206871648345162, 0.14347849887004038, -0.050911195996747316, 1.314023802253438, 0.6348258970402827, -0, 0.6109615338552246, 0.476457578966189, -2.7311218710244463, -0.052470638647269585, -0.7998764021535656, -0.7253374879437113, 0, 0.6206871648345162, -0.14347849887004038, 0.050911195996747316, -1.314023802253438, -0.6348258970402827, 0, -0.6109615338552246, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7001209024399848, -0, -0, 0.31653724289408824, -0, -0, -0.9897104202085316, 1.4350469221322282, -0.13735037357335078, -0, 0.08279503413614919, -0, -0.2797437186631715, 0.7001209024399848, 0, 0, -0.31653724289408824, 0, 0, 0.9897104202085316, -1.4350469221322282, 0.13735037357335078, 0, -0.08279503413614919, 0, 0.2797437186631715, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.7027502615621812, -0, 0.8332519091897401, -0, -0, -0.7567988418466963, -0, -0.9730343818431852, -0.5101161104860161, -0.2515400440591492, -0.16669071891829756, -0, -0.8356683570259136, 0.7027502615621812, 0, -0.8332519091897401, 0, 0, 0.7567988418466963, 0, 0.9730343818431852, 0.5101161104860161, 0.2515400440591492, 0.16669071891829756, 0, 0.8356683570259136, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, 0.7705609478497938, 0.7762132098737671, -0, -0, 0.3304567624749696, 0.24703471683023603, -0, -0, 0.2058339272568599, -0, 0.49634739711260845, 1.8970426594969672, 0, -0.7705609478497938, -0.7762132098737671, 0, 0, -0.3304567624749696, -0.24703471683023603, 0, 0, -0.2058339272568599, 0, -0.49634739711260845, -1.8970426594969672, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7354444798695736, -0, -0, 0.17797227766447388, -0, 0.8385927625193501, -0, 2.326429855201535, -0, -0, -0.6208413721884664, -0, -0.4095711945594497, -0.7354444798695736, 0, 0, -0.17797227766447388, 0, -0.8385927625193501, 0, -2.326429855201535, 0, 0, 0.6208413721884664, 0, 0.4095711945594497, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, 1.2702427184954932, 0.4206278126213235, -0.0021606686026376387, -2.382860345431941, -0, -1.4790634590370297, -0, -0, -0, 0.45831284820623924, 0, 0, 0, -1.2702427184954932, -0.4206278126213235, 0.0021606686026376387, 2.382860345431941, 0, 1.4790634590370297, 0, 0, 0, -0.45831284820623924, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.476457578966189, -2.7311218710244463, -0.052470638647269585, -0.7998764021535656, -0.7253374879437113, -0, 0.6206871648345162, -0.14347849887004038, 0.050911195996747316, -1.314023802253438, -0.6348258970402827, -0, -0.6109615338552246, -0.476457578966189, 2.7311218710244463, 0.052470638647269585, 0.7998764021535656, 0.7253374879437113, 0, -0.6206871648345162, 0.14347849887004038, -0.050911195996747316, 1.314023802253438, 0.6348258970402827, 0, 0.6109615338552246, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -1, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1.8446087238391438, 1, -0.9898577208292023, 1, 1, 2.230107839590404, 1, 0.7943912888763849, 1, 1, 1, 1, 1, -1.8446087238391438, -1, 0.9898577208292023, -1, -1, -2.230107839590404, -1, -0.7943912888763849, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}),
+			b: []float64{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},
+			c: []float64{-0.8446087238391436, 0, 1.9898577208292023, 0, 0, -2.230107839590404, 0, 0.20560871112361512, 0, 0, 0, 0, 0, 0.8446087238391436, -0, -1.9898577208292023, -0, -0, 2.230107839590404, -0, -0.20560871112361512, -0, -0, -0, -0, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+		},
+		{
+			// Dense case that panicked.
+			A: mat64.NewDense(6, 15, []float64{0.3279477313560112, 0.04126296122557327, 0.24121743535067522, -0.8676933623438741, -0.3279477313560112, -0.04126296122557327, -0.24121743535067522, 0.8676933623438741, 1, 0, 0, 0, 0, 0, 1.3702148909442915, 0.43713186538468607, 0.8613818492485417, -0.9298615442657688, -0.037784779008231184, -0.43713186538468607, -0.8613818492485417, 0.9298615442657688, 0.037784779008231184, 0, 1, 0, 0, 0, 0, 0.3478112701177931, -0, 0.748352668598051, -0.4294796840343912, -0, 0, -0.748352668598051, 0.4294796840343912, 0, 0, 0, 1, 0, 0, 0, -1, -0, 1.1913912184457485, 1.732132186658447, 0.4026384828544584, 0, -1.1913912184457485, -1.732132186658447, -0.4026384828544584, 0, 0, 0, 1, 0, 0, -0.4598555419763902, -0, 1.2088959976921831, -0.7297794575275871, 1.9835614149566971, 0, -1.2088959976921831, 0.7297794575275871, -1.9835614149566971, 0, 0, 0, 0, 1, 0, 0.5986809560819324, 0.19738159369304414, -1.0647198575836367, -0, -0.7264943883762761, -0.19738159369304414, 1.0647198575836367, 0, 0.7264943883762761, 0, 0, 0, 0, 0, 1, 0.36269644970561576}),
+			b: []float64{2.3702148909442915, 1.3478112701177931, 0, -0.4598555419763902, 1.5986809560819324, 1.3626964497056158},
+			c: []float64{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},
+		},
+		{
+			A:            mat64.NewDense(6, 11, []float64{-0.036551083288733854, -0.8967234664797694, 0.036551083288733854, 0.8967234664797694, 1, 0, 0, 0, 0, 0, -0.719908817815329, -1.9043311904524263, -0, 1.9043311904524263, 0, 0, 1, 0, 0, 0, 0, -1.142213296802784, -0, 0.17584914855696687, 0, -0.17584914855696687, 0, 0, 1, 0, 0, 0, -0.5423586338987796, -0.21663357118058713, -0.4815354890024489, 0.21663357118058713, 0.4815354890024489, 0, 0, 0, 1, 0, 0, -0.6864090947259134, -0, -0, 0, 0, 0, 0, 0, 0, 1, 0, 0.4091621839837596, -1.1853040616164046, -0.11374085137543871, 1.1853040616164046, 0.11374085137543871, 0, 0, 0, 0, 0, 1, -1.7416078575675549}),
+			b:            []float64{0.28009118218467105, -0.14221329680278405, 0.4576413661012204, 0.3135909052740866, 1.4091621839837596, -1.7416078575675549},
+			c:            []float64{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},
+			initialBasic: []int{10, 8, 7, 6, 5, 4},
+		},
+		{
+			A:            mat64.NewDense(6, 10, []float64{-0.036551083288733854, -0.8967234664797694, 0.036551083288733854, 0.8967234664797694, 1, 0, 0, 0, 0, 0, -1.9043311904524263, -0, 1.9043311904524263, 0, 0, 1, 0, 0, 0, 0, -0, 0.17584914855696687, 0, -0.17584914855696687, 0, 0, 1, 0, 0, 0, -0.21663357118058713, -0.4815354890024489, 0.21663357118058713, 0.4815354890024489, 0, 0, 0, 1, 0, 0, -0, -0, 0, 0, 0, 0, 0, 0, 1, 0, -1.1853040616164046, -0.11374085137543871, 1.1853040616164046, 0.11374085137543871, 0, 0, 0, 0, 0, 1}),
+			b:            []float64{0.28009118218467105, -0.14221329680278405, 0.4576413661012204, 0.3135909052740866, 1.4091621839837596, -1.7416078575675549},
+			c:            []float64{-1.1951160054922971, -1.354633418345746, 1.1951160054922971, 1.354633418345746, 0, 0, 0, 0, 0, 0},
+			initialBasic: []int{0, 8, 7, 6, 5, 4},
+		},
+		{
+			A:            mat64.NewDense(6, 14, []float64{-0.4398035705048233, -0, -1.1190414559968929, -0, 0.4398035705048233, 0, 1.1190414559968929, 0, 1, 0, 0, 0, 0, 0, -0, 0.45892918156139395, -0, -0, 0, -0.45892918156139395, 0, 0, 0, 1, 0, 0, 0, 0, -0, -0, -0.3163051515958635, -0, 0, 0, 0.3163051515958635, 0, 0, 0, 1, 0, 0, 0, -0, -0, -1.8226051692445888, -0.8154477101733032, 0, 0, 1.8226051692445888, 0.8154477101733032, 0, 0, 0, 1, 0, 0, -0, 1.0020104354806922, -2.80863692523519, -0.8493721031516384, 0, -1.0020104354806922, 2.80863692523519, 0.8493721031516384, 0, 0, 0, 0, 1, 0, -0.8292937871394104, -1.4615144665021647, -0, -0, 0.8292937871394104, 1.4615144665021647, 0, 0, 0, 0, 0, 0, 0, 1}),
+			b:            []float64{0, -1.0154749704172474, 0, 0, 0, -1.5002324315812783},
+			c:            []float64{1.0665389045026794, 0.097366273706136, 0, 2.7928153636989954, -1.0665389045026794, -0.097366273706136, -0, -2.7928153636989954, 0, 0, 0, 0, 0, 0},
+			initialBasic: []int{5, 12, 11, 10, 0, 8},
+		},
+		{
+			// Bad Phase I setup.
+			A: mat64.NewDense(6, 7, []float64{1.4009742075419371, 0, 0.05737255493210325, -2.5954004393412915, 0, 1.561789236911904, 0, 0.17152506517602673, 0, 0, 0, 0, 0, -0.3458126550149948, 1.900744052464951, -0.32773164134097343, -0.9648201331251137, 0, 0, 0, 0, -1.3229549190526497, 0.0692227703722903, 0, 0, -0.1024297720479933, 0.4550740188869777, 0, 0.013599438965679167, 0, 0, 0, 0, 0, -0.1164365105021209, 0, 0, 0.4077091957443405, 1.5682816151954875, 0.8411734682369051, 0.22379142247562167, 1.2206581060250778}),
+			b: []float64{0.3293809220378252, 0, -0.5688424847664554, 0, 0, 1.4832526082339592},
+			c: []float64{0.5246370956983506, -0.36608819899109946, 1.5854141981237713, 0.5170486527020665, 0, 1.4006819866163691, 0.7733814538809437},
+		},
+		{
+			// The problem is feasible, but the PhaseI problem keeps the last
+			// variable in the basis.
+			A:            mat64.NewDense(2, 3, []float64{0.7171320440380402, 0, 0.22818288617480836, 0, -0.10030202006494793, -0.3282372661549324}),
+			b:            []float64{0.8913013436978257, 0},
+			c:            []float64{0, 0, 1.16796158316812},
+			initialBasic: nil,
+		},
+		{
+			// Case where primal was returned as feasible, but the dual was returned
+			// as infeasible. This is the dual.
+			// Here, the phase I problem returns the value in the basis but equal
+			// to epsilon and not 0.
+			A: mat64.NewDense(5, 11, []float64{0.48619717875196006, 0.5089083769874058, 1.4064796473022745, -0.48619717875196006, -0.5089083769874058, -1.4064796473022745, 1, 0, 0, 0, 0, 1.5169837857318682, -0, -0, -1.5169837857318682, 0, 0, 0, 1, 0, 0, 0, -1.3096160896447528, 0.12600426735917414, 0.296082394213142, 1.3096160896447528, -0.12600426735917414, -0.296082394213142, 0, 0, 1, 0, 0, -0, -0, 1.9870800277141467, 0, 0, -1.9870800277141467, 0, 0, 0, 1, 0, -0.3822356988571877, -0, -0.1793908926957139, 0.3822356988571877, 0, 0.1793908926957139, 0, 0, 0, 0, 1}),
+			b: []float64{0.6015865977347667, 0, -1.5648780993757594, 0, 0},
+			c: []float64{-0.642801659201449, -0.5412741400343285, -1.4634460998530177, 0.642801659201449, 0.5412741400343285, 1.4634460998530177, 0, 0, 0, 0, 0},
+		},
+		{
+			// Caused linear solve error. The error is because replacing the minimum
+			// index in Bland causes the new basis to be singular. This
+			// necessitates the ending loop in bland over possible moves.
+			A: mat64.NewDense(9, 23, []float64{-0.898219823758102, -0, -0, -0, 1.067555075209233, 1.581598470243863, -1.0656096883610071, 0.898219823758102, 0, 0, 0, -1.067555075209233, -1.581598470243863, 1.0656096883610071, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1.5657353278668433, 0.5798888118401012, -0, 0.14560553520321928, -0, -0, -0, 1.5657353278668433, -0.5798888118401012, 0, -0.14560553520321928, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -0, -0, -1.5572250142582087, -0, -0, -0, -0, 0, 0, 1.5572250142582087, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -0, -0, -0, -1.1266215512973428, -0, 1.0661059397023553, -0, 0, 0, 0, 1.1266215512973428, 0, -1.0661059397023553, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -0, -2.060232129551813, 1.756900609902372, -0, -0, -0, -0, 0, 2.060232129551813, -1.756900609902372, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1.0628806512935949, -0, -0, 0.3306985942820342, -0, 0.5013194822231914, -0, -1.0628806512935949, 0, 0, -0.3306985942820342, 0, -0.5013194822231914, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -0.02053916418367785, 2.0967009672108627, -0, 1.276296057052031, -0, -0.8396554873675388, -0, 0.02053916418367785, -2.0967009672108627, 0, -1.276296057052031, 0, 0.8396554873675388, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1.5173172721095745, -0, -0, -0, -0, -0.7781977786718928, -0.08927683907374018, 1.5173172721095745, 0, 0, 0, 0, 0.7781977786718928, 0.08927683907374018, 0, 0, 0, 0, 0, 0, 0, 1, 0, -0, -0, -0, -0, -0, 0.39773149008355624, -0, 0, 0, 0, 0, 0, -0.39773149008355624, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}),
+			b: []float64{0, 0, 0, 0, 0, 0, 0, 0, 0},
+			c: []float64{0.24547850255842107, -0.9373919913433648, 0, 0, 0, 0.2961224049153204, 0, -0.24547850255842107, 0.9373919913433648, -0, -0, -0, -0.2961224049153204, -0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+		},
+		{
+			// Caused error because ALL of the possible replacements in Bland cause.
+			// ab to be singular. This necessitates outer loop in bland over possible
+			// moves.
+			A:            mat64.NewDense(9, 23, []float64{0.6595219196440785, -0, -0, -1.8259394918781682, -0, -0, 0.005457361044175046, -0.6595219196440785, 0, 0, 1.8259394918781682, 0, 0, -0.005457361044175046, 1, 0, 0, 0, 0, 0, 0, 0, 0, -0, -0.10352878714214864, -0, -0, -0, -0, 0.5945016966696087, 0, 0.10352878714214864, 0, 0, 0, 0, -0.5945016966696087, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.31734882842876444, -0, -0, -0, -0, -0, -0.716633126367685, -0.31734882842876444, 0, 0, 0, 0, 0, 0.716633126367685, 0, 0, 1, 0, 0, 0, 0, 0, 0, -0.7769812182932578, -0, -0.17370050158829553, 0.19405062263734607, -0, 1.1472330031002533, -0.6776631768730962, 0.7769812182932578, 0, 0.17370050158829553, -0.19405062263734607, 0, -1.1472330031002533, 0.6776631768730962, 0, 0, 0, 1, 0, 0, 0, 0, 0, -0, 0.8285611486611473, -0, -0, -0, -0, -0, 0, -0.8285611486611473, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -2.088953453647358, 1.3286488791152795, -0, -0, -0, -0, 0.9147833235021142, 2.088953453647358, -1.3286488791152795, 0, 0, 0, 0, -0.9147833235021142, 0, 0, 0, 0, 0, 1, 0, 0, 0, -0, -0, -0, -0, 0.6560365621262937, -0, -0, 0, 0, 0, 0, -0.6560365621262937, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.8957188338098074, -0, -0, -0, -0, -0, -0, -0.8957188338098074, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, -0.2761381891117365, -0, -0, -0, 1.1154921426237823, 0.06429872020552618, -0, 0.2761381891117365, 0, 0, 0, -1.1154921426237823, -0.06429872020552618, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}),
+			b:            []float64{0, 0, 0, 0, 0.5046208538522362, 1.0859412982429362, -2.066283584195025, 0, -0.2604305274353169},
+			c:            []float64{0, 0, 0, 0, 0, 0, -0.05793762969330718, -0, -0, -0, -0, -0, -0, 0.05793762969330718, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+			initialBasic: []int{22, 11, 7, 19, 18, 17, 16, 15, 14},
+		},
+		{
+			// Caused initial supplied basis of Phase I to be singular.
+			A: mat64.NewDense(7, 11, []float64{0, 0, 0, 0, 0, 0.6667874223914787, -0.04779440888372957, -0.810020924434026, 0, 1.4190243477163373, 0, 0, 1.0452496826112936, 1.1966134226828076, 0, 0, 0, 0, -0.676136041089015, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.123232807871834, 0.2795467733707712, 0.21997115467272987, 0, -0.1572003980840453, 0, 0, 0, 0, 0.5130196002804861, 0, -0.005957174211761673, 0.3262874931735277, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.5582052881594286, 0, 0.3544026193217651, 0, -1.0761986709145068, 0, 0.2438593072108347, 0, 0, 0, 0, 1.387509848081664, 0, 0, 0.3958750570508226, 1.6281679612990678, 0, 0, -0.24638311667922103, 0, 0, 0, 0, 0, 0, -0.628850893994423}),
+			b: []float64{0.4135281763629115, 0, 0, 0, 0, 0, 0},
+			c: []float64{0.5586772876113472, 0, 0.14261332948424457, 0, -0.016394076753000086, -0.506087285562544, 0, 0.37619482505459145, 1.2943822852419233, 0.5887960293578207, 0},
+		},
+	} {
+		testSimplex(t, test.initialBasic, test.c, test.A, test.b, convergenceTol)
+	}
+
+	rnd := rand.New(rand.NewSource(1))
+	// Randomized tests
+	testRandomSimplex(t, 20000, 0.7, 10, rnd)
+	testRandomSimplex(t, 20000, 0, 10, rnd)
+	testRandomSimplex(t, 200, 0, 100, rnd)
+	testRandomSimplex(t, 2, 0, 400, rnd)
+}
+
+func testRandomSimplex(t *testing.T, nTest int, pZero float64, maxN int, rnd *rand.Rand) {
+	// Try a bunch of random LPs
+	for i := 0; i < nTest; i++ {
+		n := rnd.Intn(maxN) + 2 // n must be at least two.
+		m := rnd.Intn(n-1) + 1  // m must be between 1 and n
+		if m == 0 || n == 0 {
+			continue
+		}
+		randValue := func() float64 {
+			//var pZero float64
+			v := rnd.Float64()
+			if v < pZero {
+				return 0
+			}
+			return rnd.NormFloat64()
+		}
+		a := mat64.NewDense(m, n, nil)
+		for i := 0; i < m; i++ {
+			for j := 0; j < n; j++ {
+				a.Set(i, j, randValue())
+			}
+		}
+		b := make([]float64, m)
+		for i := range b {
+			b[i] = randValue()
+		}
+
+		c := make([]float64, n)
+		for i := range c {
+			c[i] = randValue()
+		}
+
+		testSimplex(t, nil, c, a, b, convergenceTol)
+	}
+}
+
+func testSimplex(t *testing.T, initialBasic []int, c []float64, a mat64.Matrix, b []float64, convergenceTol float64) error {
+	primalOpt, primalX, _, errPrimal := simplex(initialBasic, c, a, b, convergenceTol)
+	if errPrimal == nil {
+		// No error solving the simplex, check that the solution is feasible.
+		var bCheck mat64.Vector
+		bCheck.MulVec(a, mat64.NewVector(len(primalX), primalX))
+		if !mat64.EqualApprox(&bCheck, mat64.NewVector(len(b), b), 1e-10) {
+			t.Errorf("No error in primal but solution infeasible")
+		}
+	}
+
+	primalInfeasible := errPrimal == ErrInfeasible
+	primalUnbounded := errPrimal == ErrUnbounded
+	primalBounded := errPrimal == nil
+	primalASingular := errPrimal == ErrSingular
+	primalZeroRow := errPrimal == ErrZeroRow
+	primalZeroCol := errPrimal == ErrZeroColumn
+
+	primalBad := !primalInfeasible && !primalUnbounded && !primalBounded && !primalASingular && !primalZeroRow && !primalZeroCol
+
+	// It's an error if it's not one of the known returned errors. If it's
+	// singular the problem is undefined and so the result cannot be compared
+	// to the dual.
+	if errPrimal == ErrSingular || primalBad {
+		if primalBad {
+			t.Errorf("non-known error returned: %s", errPrimal)
+		}
+		return errPrimal
+	}
+
+	// Compare the result to the answer found from solving the dual LP.
+
+	// Construct and solve the dual LP.
+	// Standard Form:
+	//  minimize c^T * x
+	//    subject to  A * x = b, x >= 0
+	// The dual of this problem is
+	//  maximize -b^T * nu
+	//   subject to A^T * nu + c >= 0
+	// Which is
+	//   minimize b^T * nu
+	//   subject to -A^T * nu <= c
+
+	negAT := &mat64.Dense{}
+	negAT.Clone(a.T())
+	negAT.Scale(-1, negAT)
+	cNew, aNew, bNew := Convert(b, negAT, c, nil, nil)
+
+	dualOpt, dualX, _, errDual := simplex(nil, cNew, aNew, bNew, convergenceTol)
+	if errDual == nil {
+		// Check that the dual is feasible
+		var bCheck mat64.Vector
+		bCheck.MulVec(aNew, mat64.NewVector(len(dualX), dualX))
+		if !mat64.EqualApprox(&bCheck, mat64.NewVector(len(bNew), bNew), 1e-10) {
+			t.Errorf("No error in dual but solution infeasible")
+		}
+	}
+
+	// Check about the zero status.
+	if errPrimal == ErrZeroRow || errPrimal == ErrZeroColumn {
+		return errPrimal
+	}
+
+	// If the primal problem is feasible, then the primal and the dual should
+	// be the same answer. We have flopped the sign in the dual (minimizing
+	// b^T *nu instead of maximizing -b^T*nu), so flip it back.
+	if errPrimal == nil {
+		if errDual != nil {
+			fmt.Println("errDual", errDual)
+			panic("here")
+			t.Errorf("Primal feasible but dual errored: %s", errDual)
+		}
+		dualOpt *= -1
+		if !floats.EqualWithinAbsOrRel(dualOpt, primalOpt, convergenceTol, convergenceTol) {
+			t.Errorf("Primal and dual value mismatch. Primal %v, dual %v.", primalOpt, dualOpt)
+		}
+	}
+	// If the primal problem is unbounded, then the dual should be infeasible.
+	if errPrimal == ErrUnbounded && errDual != ErrInfeasible {
+		t.Errorf("Primal unbounded but dual not infeasible. ErrDual = %s", errDual)
+	}
+
+	// If the dual is unbounded, then the primal should be infeasible.
+	if errDual == ErrUnbounded && errPrimal != ErrInfeasible {
+		t.Errorf("Dual unbounded but primal not infeasible. ErrDual = %s", errPrimal)
+	}
+
+	// If the primal is infeasible, then the dual should be either infeasible
+	// or unbounded.
+	if errPrimal == ErrInfeasible {
+		if errDual != ErrUnbounded && errDual != ErrInfeasible && errDual != ErrZeroColumn {
+			t.Errorf("Primal infeasible but dual not infeasible or unbounded: %s", errDual)
+		}
+	}
+
+	return errPrimal
+}
diff --git a/convex/lp/simplexexample_test.go b/convex/lp/simplexexample_test.go
new file mode 100644
index 0000000..b71fb7e
--- /dev/null
+++ b/convex/lp/simplexexample_test.go
@@ -0,0 +1,29 @@
+// 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 lp_test
+
+import (
+	"fmt"
+	"log"
+
+	"github.com/gonum/matrix/mat64"
+	"github.com/gonum/optimize/convex/lp"
+)
+
+func Example_Simplex() {
+	c := []float64{-1, -2, 0, 0}
+	A := mat64.NewDense(2, 4, []float64{-1, 2, 1, 0, 3, 1, 0, 1})
+	b := []float64{4, 9}
+
+	opt, x, err := lp.Simplex(c, A, b, 0, nil)
+	if err != nil {
+		log.Fatal(err)
+	}
+	fmt.Printf("opt: %v\n", opt)
+	fmt.Printf("x: %v\n", x)
+	// Output:
+	// opt: -8
+	// x: [2 3 0 0]
+}