diff --git a/Control/Monad/LPMonad/Internal.hs b/Control/Monad/LPMonad/Internal.hs
--- a/Control/Monad/LPMonad/Internal.hs
+++ b/Control/Monad/LPMonad/Internal.hs
@@ -40,14 +40,16 @@
 	varBds,
 	setVarBounds,
 	setVarKind,
-	newVariables,
-	newVariables') where
+-- 	newVariables,
+-- 	newVariables'
+	) where
 
 import Control.Monad.State.Strict
 import Control.Monad.Identity
 
 import Data.Map
 import Data.Monoid
+import Data.Functor
 
 import Data.LinearProgram.Common
 
@@ -91,7 +93,11 @@
 	(Ord v, Group c, Monad m) => LinFunc v c -> LinFunc v c -> LPT v c m () #-}
 {-# SPECIALIZE geq :: (Ord v, Group c) => LinFunc v c -> LinFunc v c -> LPM v c (),
 	(Ord v, Group c, Monad m) => LinFunc v c -> LinFunc v c -> LPT v c m () #-}
--- | Specifies the relationship between two functions in the variables.
+-- | Specifies the relationship between two functions in the variables.  So, for example,
+-- 
+-- > equal (f ^+^ g) h
+-- 
+-- constrains the value of @h@ to be equal to the value of @f@ plus the value of @g@.
 equal, leq, geq :: (Ord v, Group c, MonadState (LP v c) m) => LinFunc v c -> LinFunc v c -> m ()
 equal f g = equalTo (f ^-^ g) zero
 leq f g = leqTo (f ^-^ g) zero
@@ -130,39 +136,39 @@
 leqTo' lab f v = constrain' lab f (UBound v)
 geqTo' lab f v = constrain' lab f (LBound v)
 
-{-# SPECIALIZE newVariables :: (Ord v, Enum v) => Int -> LPM v c [v],
-	(Ord v, Enum v, Monad m) => Int -> LPT v c m [v] #-}
--- | Returns a list of @k@ unused variables.  If the program is currently empty,
--- starts at @'toEnum' 0@.  Otherwise, if @v@ is the biggest variable currently in use
--- (by the 'Ord' ordering), then this returns @take k (tail [v..])@, which uses the 'Enum'
--- implementation.  Note that if the 'Enum' instance doesn't play well with 'Ord',
--- bad things can happen.
-newVariables :: (MonadState (LP v c) m, Ord v, Enum v) => Int -> m [v]
-newVariables !k = do	LP{..} <- get
-			let allVars0 = fmap (const ()) objective `union`
-				unions [fmap (const ()) f | Constr _ f _ <- constraints] `union`
-				fmap (const ()) varBounds `union` fmap (const ()) varTypes
-			case minViewWithKey allVars0 of
-				Nothing	-> return $ take k [toEnum 0..]
-				Just ((start, _), _)
-					-> return $ take k $ tail [start..]
-					
-{-# SPECIALIZE newVariables' :: (Ord v, Enum v) => LPM v c [v],
-	(Ord v, Enum v, Monad m) => LPT v c m [v] #-}
--- | Returns an infinite list of unused variables.  If the program is currently empty,
--- starts at @'toEnum' 0@.  Otherwise, if @v@ is the biggest variable currently in use
--- (by the 'Ord' ordering), then this returns @tail [v..]@, which uses the 'Enum'
--- implementation.  Note that if the 'Enum' instance doesn't play well with 'Ord',
--- bad things can happen.
-newVariables' :: (MonadState (LP v c) m, Ord v, Enum v) => m [v]
-newVariables' = do	LP{..} <- get
-			let allVars0 = fmap (const ()) objective `union`
-				unions [fmap (const ()) f | Constr _ f _ <- constraints] `union`
-				fmap (const ()) varBounds `union` fmap (const ()) varTypes
-			case minViewWithKey allVars0 of
-				Nothing	-> return [toEnum 0..]
-				Just ((start, _), _)
-					-> return $ tail [start..]
+-- {-# SPECIALIZE newVariables :: (Ord v, Enum v) => Int -> LPM v c [v],
+-- 	(Ord v, Enum v, Monad m) => Int -> LPT v c m [v] #-}
+-- -- | Returns a list of @k@ unused variables.  If the program is currently empty,
+-- -- starts at @'toEnum' 0@.  Otherwise, if @v@ is the biggest variable currently in use
+-- -- (by the 'Ord' ordering), then this returns @take k (tail [v..])@, which uses the 'Enum'
+-- -- implementation.  Note that if the 'Enum' instance doesn't play well with 'Ord',
+-- -- bad things can happen.
+-- newVariables :: (MonadState (LP v c) m, Ord v, Enum v) => Int -> m [v]
+-- newVariables !k = do	LP{..} <- get
+-- 			let allVars0 = () <$ objective `union`
+-- 				unions [() <$ f | Constr _ f _ <- constraints] `union`
+-- 				(() <$ varBounds) `union` (() <$ varTypes)
+-- 			case minViewWithKey allVars0 of
+-- 				Nothing	-> return $ take k [toEnum 0..]
+-- 				Just ((start, _), _)
+-- 					-> return $ take k $ tail [start..]
+-- 					
+-- {-# SPECIALIZE newVariables' :: (Ord v, Enum v) => LPM v c [v],
+-- 	(Ord v, Enum v, Monad m) => LPT v c m [v] #-}
+-- -- | Returns an infinite list of unused variables.  If the program is currently empty,
+-- -- starts at @'toEnum' 0@.  Otherwise, if @v@ is the biggest variable currently in use
+-- -- (by the 'Ord' ordering), then this returns @tail [v..]@, which uses the 'Enum'
+-- -- implementation.  Note that if the 'Enum' instance doesn't play well with 'Ord',
+-- -- bad things can happen.
+-- newVariables' :: (MonadState (LP v c) m, Ord v, Enum v) => m [v]
+-- newVariables' = do	LP{..} <- get
+-- 			let allVars0 = () <$ objective `union`
+-- 				unions [() <$ f | Constr _ f _ <- constraints] `union`
+-- 				(() <$ varBounds) `union` (() <$ varTypes)
+-- 			case minViewWithKey allVars0 of
+-- 				Nothing	-> return [toEnum 0..]
+-- 				Just ((start, _), _)
+-- 					-> return $ tail [start..]
 
 {-# SPECIALIZE varEq :: (Ord v, Ord c) => v -> c -> LPM v c (),
 	(Ord v, Ord c, Monad m) => v -> c -> LPT v c m () #-}
@@ -215,7 +221,7 @@
 addObjective :: (Ord v, Group c, MonadState (LP v c) m) => LinFunc v c -> m ()
 addObjective obj = modify addObj where
 	addObj lp@LP{..} = lp {objective = obj ^+^ objective}
-		
+
 {-# SPECIALIZE addWeightedObjective :: (Ord v, Module r c) => r -> LinFunc v c -> LPM v c (),
 	(Ord v, Module r c, Monad m) => r -> LinFunc v c -> LPT v c m () #-}
 -- | Adds this function to the objective function, with the specified weight.  Equivalent to
diff --git a/Data/Algebra.hs b/Data/Algebra.hs
--- a/Data/Algebra.hs
+++ b/Data/Algebra.hs
@@ -12,7 +12,9 @@
 	-- * Algebraic structures
 	Group(..),
 	Ring(..),
+	Field(..),
 	Module(..),
+	VectorSpace(..),
 	Poly,
 	varPoly,
 	GroupRing,
@@ -29,4 +31,5 @@
 
 import Data.Algebra.Group
 import Data.Algebra.Ring
+import Data.Algebra.Field
 import Data.Algebra.Module
diff --git a/Data/Algebra/Field.hs b/Data/Algebra/Field.hs
new file mode 100644
--- /dev/null
+++ b/Data/Algebra/Field.hs
@@ -0,0 +1,24 @@
+{-# LANGUAGE UndecidableInstances, FlexibleInstances, MultiParamTypeClasses #-}
+
+module Data.Algebra.Field where
+
+import Data.Ratio
+
+import Data.Algebra.Ring
+import Data.Algebra.Module
+
+class Ring f => Field f where
+	inv :: f -> f
+	(/#) :: f -> f -> f
+	inv x = one /# x
+	a /# b = a *# inv b
+
+instance Field Double where
+	inv = recip
+
+instance Integral a => Field (Ratio a) where
+	{-# SPECIALIZE instance Field Rational #-}
+	inv = recip
+
+class (Module f v, Field f) => VectorSpace f v
+instance (Module f v, Field f) => VectorSpace f v
diff --git a/Data/Algebra/Module.hs b/Data/Algebra/Module.hs
--- a/Data/Algebra/Module.hs
+++ b/Data/Algebra/Module.hs
@@ -76,11 +76,11 @@
 type LinFunc = M.Map
 
 -- | Given a variable @v@, returns the function equivalent to @v@.
-var :: (Ord v, Num c) => v -> LinFunc v c
-var v = M.singleton v 1
+var :: (Ord v, Ring c) => v -> LinFunc v c
+var v = M.singleton v one
 
 -- | @c '*&' v@ is equivalent to @c '*^' 'var' v@.
-(*&) :: (Ord v, Num c) => c -> v -> LinFunc v c
+(*&) :: (Ord v, Ring c) => c -> v -> LinFunc v c
 c *& v = M.singleton v c
 
 -- | Equivalent to @'vsum' . 'map' 'var'@.
diff --git a/Data/Algebra/Ring.hs b/Data/Algebra/Ring.hs
--- a/Data/Algebra/Ring.hs
+++ b/Data/Algebra/Ring.hs
@@ -40,15 +40,18 @@
 	one = 1
 	(*#) = (*)
 
+-- | The polynomial ring.
 instance Ring r => Ring (Poly r) where
 	one = [one]
 	(p:ps) *# (q:qs) = (p *# q):(ps *# (q:qs) ^+^ map (p *#) qs)
 	_ *# _ = []
 
+-- | The function ring.
 instance Ring r => Ring (a -> r) where
 	one = const one
 	(*#) = liftA2 (*#)
 
+-- | The group ring.
 instance (Ord g, Group g, Ring r) => Ring (GroupRing r g) where
 	one = M.singleton zero one
 	m *# n = M.fromListWith (^+^) [(u ^+^ v, f *# g) | (u, f) <- M.assocs m, (v, g) <- M.assocs n]
diff --git a/Data/LinearProgram/GLPK/IO/Internal.hs b/Data/LinearProgram/GLPK/IO/Internal.hs
--- a/Data/LinearProgram/GLPK/IO/Internal.hs
+++ b/Data/LinearProgram/GLPK/IO/Internal.hs
@@ -6,6 +6,8 @@
 import Control.Monad.Trans (liftIO, lift)
 
 import Data.Map hiding (map, filter)
+import Debug.Trace
+import Foreign.Storable
 
 import Data.LinearProgram.Common
 import Data.LinearProgram.GLPK.Common
@@ -78,7 +80,7 @@
 	case typ of
 		1	-> return Free
 		2	-> liftM LBound (lb i)
-		3	-> liftM UBound (lb i)
+		3	-> liftM UBound (ub i)
 		4	-> liftM2 Bound (lb i) (ub i)
 		_	-> liftM Equ (lb i)
 		
@@ -88,12 +90,12 @@
 getRows :: GLPK [(Int, [(Int, Double)])]
 getRows = do	n <- getNumRows
 		m <- getNumCols
-		ixs <- liftIO $ mallocArray m
-		coefs <- liftIO $ mallocArray m
+		ixs <- liftIO $ mallocArray (m+1)
+		coefs <- liftIO $ mallocArray (m+1)
 		sequence [do
 			k <- liftM fromIntegral $ GLP $ \ lp -> glpGetMatRow lp (fromIntegral i) ixs coefs
-			ixsL <- liftIO $ peekArray k ixs
-			coefsL <- liftIO $ peekArray k coefs
+			ixsL <- liftIO $ mapM (peekElemOff ixs) [1..k]
+			coefsL <- liftIO $ mapM (peekElemOff ixs) [1..k]
 			return (i, zip (map fromIntegral ixsL) (map realToFrac coefsL))
 			| i <- [1..n]]
 
diff --git a/Data/LinearProgram/LinExpr.hs b/Data/LinearProgram/LinExpr.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/LinExpr.hs
@@ -0,0 +1,58 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+module Data.LinearProgram.LinExpr (LinExpr(..), solve, substituteExpr, simplifyExpr,
+	constTerm, coeffTerm, funcToExpr) where
+
+import Control.Monad
+
+import Data.LinearProgram.Types
+import Data.Algebra
+import Data.Functor
+import Data.Foldable
+
+import Data.Map
+
+import Prelude hiding (lookup, filter, foldr)
+
+constTerm :: LinExpr v c -> c
+constTerm (LinExpr _ c) = c
+
+coeffTerm :: LinExpr v c -> LinFunc v c
+coeffTerm (LinExpr a _) = a
+
+funcToExpr :: Group c => LinFunc v c -> LinExpr v c
+funcToExpr = flip LinExpr zero
+
+data LinExpr v c = LinExpr (LinFunc v c) c deriving (Eq, Read, Show)
+
+instance (Ord v, Group c) => Group (LinExpr v c) where
+	zero = LinExpr zero zero
+	LinExpr a1 c1 ^+^ LinExpr a2 c2 = LinExpr (a1 ^+^ a2) (c1 ^+^ c2)
+	LinExpr a1 c1 ^-^ LinExpr a2 c2 = LinExpr (a1 ^-^ a2) (c1 ^-^ c2)
+	neg (LinExpr a c) = LinExpr (neg a) (neg c)
+
+instance (Ord v, Module r c) => Module r (LinExpr v c) where
+	k *^ LinExpr a c = LinExpr (k *^ a) (k *^ c)
+
+substituteExpr :: (Ord v, Module c c) => v -> LinExpr v c -> LinExpr v c -> LinExpr v c
+substituteExpr v expV expr@(LinExpr a c) = case lookup v a of
+	Nothing	-> expr
+	Just k	-> LinExpr (delete v a) c ^+^ (k *^ expV)
+
+simplifyExpr :: (Ord v, Module c c) => LinExpr v c -> Map v (LinExpr v c) -> LinExpr v c
+simplifyExpr (LinExpr a c) sol =
+	foldrWithKey (const (^+^)) (LinExpr (difference a sol) c) (intersectionWith (*^) a sol)
+
+solve :: (Ord v, Eq c, VectorSpace c c) => [(LinFunc v c, c)] -> Maybe (Map v (LinExpr v c))
+solve equs = solve' [LinExpr a (neg c) | (a, c) <- equs]
+
+solve' :: (Ord v, Eq c, VectorSpace c c) => [LinExpr v c] -> Maybe (Map v (LinExpr v c))
+solve' (LinExpr a c:equs) = case minViewWithKey (filter (/= zero) a) of
+	Nothing	-> guard (c == zero) >> solve' equs
+	Just ((x, a0), a') -> let expX = neg (inv a0 *^ LinExpr a' c) in
+		liftM (simplifyExpr expX >>= insert x) (solve' (substituteExpr x expX <$> equs))
+solve' [] = return empty
+
+{-# RULES
+	"mapWithKey/mapWithKey" forall f g m .
+		mapWithKey f (mapWithKey g m) = mapWithKey (liftM2 (.) f g) m
+	#-}
diff --git a/Data/LinearProgram/Spec.hs b/Data/LinearProgram/Spec.hs
--- a/Data/LinearProgram/Spec.hs
+++ b/Data/LinearProgram/Spec.hs
@@ -1,12 +1,14 @@
 {-# LANGUAGE TupleSections, RecordWildCards, DeriveFunctor #-}
 module Data.LinearProgram.Spec (Constraint(..), VarTypes, ObjectiveFunc, VarBounds, LP(..),
-	mapVars, mapVals) where
+	mapVars, mapVals, allVars) where
 
 import Control.Applicative ((<$>))
+import Control.DeepSeq
 import Control.Monad
 
 import Data.Char (isSpace)
 import Data.Monoid
+import Data.Functor
 import Data.Map hiding (map)
 
 import Text.ParserCombinators.ReadP
@@ -31,6 +33,10 @@
 data LP v c = LP {direction :: Direction, objective :: ObjectiveFunc v c, constraints :: [Constraint v c],
 			varBounds :: VarBounds v c, varTypes :: VarTypes v} deriving (Read, Show, Functor)
 
+allVars :: Ord v => LP v c -> Map v ()
+allVars LP{..} = foldl union ((() <$ objective) `union` (() <$ varBounds) `union` (() <$ varTypes))
+	[() <$ f | Constr _ f _ <- constraints]
+
 showBds :: Show c => String -> Bounds c -> String
 showBds expr bds = case bds of
 	Free	-> expr ++ " free"
@@ -140,9 +146,20 @@
 mapVals :: (c -> c') -> LP v c -> LP v c'
 mapVals = fmap
 
--- instance (NFData v, NFData c) => NFData (Constraint v c) where
--- 	rnf (Constr lab f b) = lab `deepseq` f `deepseq` rnf b
+instance (NFData v, NFData c) => NFData (Constraint v c) where
+	rnf (Constr lab f b) = lab `deepseq` f `deepseq` rnf b
 
--- instance (NFData v, NFData c) => NFData (LP v c) where
--- 	rnf LP{..} = direction `deepseq` objective `deepseq` constraints `deepseq`
--- 		varBounds `deepseq` rnf varTypes
+instance (NFData v, NFData c) => NFData (LP v c) where
+	rnf LP{..} = direction `deepseq` objective `deepseq` constraints `deepseq`
+		varBounds `deepseq` rnf varTypes
+
+instance NFData VarKind
+instance NFData Direction
+instance NFData c => NFData (Bounds c) where
+	rnf Free = ()
+	rnf (Equ c) = rnf c
+	rnf (LBound c) = rnf c
+	rnf (UBound c) = rnf c
+	rnf (Bound l u) = l `deepseq` rnf u
+-- instance (NFData k, NFData a) => NFData (Map k a) where
+-- 	rnf m = foldrWithKey (\ k a -> deepseq k . deepseq a) () m
diff --git a/Data/LinearProgram/Types.hs b/Data/LinearProgram/Types.hs
--- a/Data/LinearProgram/Types.hs
+++ b/Data/LinearProgram/Types.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE DeriveFunctor #-}
 module Data.LinearProgram.Types (VarKind(..), Direction(..), Bounds(..)) where
 
--- import Control.DeepSeq
+import Control.DeepSeq
 
 import Data.Monoid
 
@@ -18,7 +18,7 @@
 -- instance NFData Direction
 
 data Bounds a =
-	Free | LBound a | UBound a | Equ a | Bound a a deriving (Eq, Show, Read, Functor)
+	Free | LBound !a | UBound !a | Equ !a | Bound !a !a deriving (Eq, Show, Read, Functor)
 
 -- instance NFData (Bounds a)
 
diff --git a/glpk-hs.cabal b/glpk-hs.cabal
--- a/glpk-hs.cabal
+++ b/glpk-hs.cabal
@@ -1,5 +1,5 @@
 Name:           glpk-hs
-Version:        0.2.4
+Version:        0.3.1
 Author:         Louis Wasserman
 License:        BSD3
 License-file:   LICENSE
@@ -23,9 +23,10 @@
 
 extra-source-files: examples/example1.hs
 
-Build-Depends:    base >= 4 && < 5, array, containers, mtl
+Build-Depends:    base >= 4 && < 5, array, containers, mtl, deepseq
 Exposed-modules:  Data.LinearProgram,
                   Data.LinearProgram.Common,
+		  Data.LinearProgram.LinExpr,
                   Data.LinearProgram.GLPK,
                   Data.LinearProgram.GLPK.Solver,
                   Data.LinearProgram.GLPK.IO,
@@ -42,6 +43,7 @@
                   Data.LinearProgram.Types,
                   Data.Algebra.Group,
                   Data.Algebra.Ring,
-                  Data.Algebra.Module
+                  Data.Algebra.Module,
+		  Data.Algebra.Field
 c-sources:        glpk/glpk.c
 extra-libraries:  glpk
diff --git a/glpk/glpk.c b/glpk/glpk.c
--- a/glpk/glpk.c
+++ b/glpk/glpk.c
@@ -105,7 +105,7 @@
 int c_glp_mip_solve(glp_prob *lp, int msg_lev, int br_tech, int bt_tech, int pp_tech,
 		     	int fp_heur, int tm_lim, int cuts, double mip_gap, int presolve){
   	glp_iocp iocp;
-	glp_mem_limit(750);
+	glp_mem_limit(10000);
 // 	printf ("%d %d %d time\n", msg_lev, br_tech, tm_lim);
 	glp_init_iocp(&iocp);
 	iocp.msg_lev = msg_lev;
@@ -194,7 +194,7 @@
 }
 
 int c_glp_get_mat_row (glp_prob *lp, int i, int ind[], double val[]){
-  	return c_glp_get_mat_row (lp, i, ind, val);
+  	return glp_get_mat_row (lp, i, ind, val);
 }
 
 int c_glp_read_lp(glp_prob *lp, const char * fname){