diff --git a/Control/Monad/LPMonad.hs b/Control/Monad/LPMonad.hs
deleted file mode 100644
--- a/Control/Monad/LPMonad.hs
+++ /dev/null
@@ -1,98 +0,0 @@
-{-# LANGUAGE FlexibleContexts #-}
-
--- | A collection of operations that can be used to specify linear programming in a
--- simple, monadic way.  It is not too difficult to construct 'LP' values explicitly,
--- but this module may help simplify and modularize the construction of the linear program,
--- for example separating different families of constraints in the problem specification.
--- 
--- Many of these functions should be executed in either the @'LPM' v c@ or the @'LPT' v c 'IO'@ monad.
--- If you wish to generate new variables on an ad-hoc basis, rather than supplying your own variable type, use the
--- 'VSupply' or 'VSupplyT' monads in your transformer stack, as in @'LPT' 'Var' c 'VSupply'@ or
--- @'LPT' 'Var' c ('VSupplyT' 'IO')@.  To generate new variables, use 'supplyNew' or 'supplyN'.
-module Control.Monad.LPMonad (
-	module Control.Monad.LPMonad.Internal,
-	-- * Generation of new variables
-	module Control.Monad.LPMonad.Supply,
-	-- * Solvers
-	quickSolveMIP,
-	quickSolveLP,
-	glpSolve,
-	quickSolveMIP',
-	quickSolveLP',
-	glpSolve',
-	-- * File I/O
-	writeLPToFile,
-	readLPFromFile,
-	readLPFromFile') where
-
-import Control.Monad ((<=<))
-import Control.Monad.State.Class (MonadState(..))
-import Control.Monad.Trans (MonadIO (..))
-
-import Data.Map (Map)
-
-import Data.LinearProgram.Common
-import Control.Monad.LPMonad.Internal
-import Control.Monad.LPMonad.Supply
-
-import Data.LinearProgram.GLPK.Solver
-import Data.LinearProgram.GLPK.IO
-
-{-# SPECIALIZE quickSolveLP :: (Ord v, Real c) => 
-	LPT v c IO (ReturnCode, Maybe (Double, Map v Double)) #-}
-{-# SPECIALIZE quickSolveMIP :: (Ord v, Real c) => 
-	LPT v c IO (ReturnCode, Maybe (Double, Map v Double)) #-}
--- | Solves the linear program with the default settings in GLPK.  Returns the return code,
--- and if the solver was successful, the objective function value and the settings of each variable.
-quickSolveLP, quickSolveMIP :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => 
-	m (ReturnCode, Maybe (Double, Map v Double))
-quickSolveLP = glpSolve simplexDefaults
-quickSolveMIP = glpSolve mipDefaults
-
-{-# SPECIALIZE glpSolve :: (Ord v, Real c) => GLPOpts -> LPT v c IO (ReturnCode, Maybe (Double, Map v Double)) #-}
--- | Solves the linear program with the specified options in GLPK.  Returns the return code,
--- and if the solver was successful, the objective function value and the settings of each variable.
-glpSolve :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => GLPOpts -> m (ReturnCode, Maybe (Double, Map v Double))
-glpSolve opts = get >>= liftIO . glpSolveVars opts
-
-{-# SPECIALIZE quickSolveLP' :: (Ord v, Real c) => LPT v c IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c])) #-}
-{-# SPECIALIZE quickSolveMIP' :: (Ord v, Real c) => LPT v c IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c])) #-}
--- | Solves the linear program with the default settings in GLPK.  Returns the return code,
--- and if the solver was successful, the objective function value, the settings of each variable, and the
--- value of each constraint/row.
-quickSolveLP', quickSolveMIP' :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => 
-	m (ReturnCode, Maybe (Double, Map v Double, [RowValue v c]))
-quickSolveLP' = glpSolve' simplexDefaults
-quickSolveMIP' = glpSolve' mipDefaults
-
-{-# SPECIALIZE glpSolve' :: (Ord v, Real c) => GLPOpts -> LPT v c IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c])) #-}
--- | Solves the linear program with the specified options in GLPK.  Returns the return code,
--- and if the solver was successful, the objective function value, the settings of each variable, and
--- the value of each constraint/row.
-glpSolve' :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => 
-	GLPOpts -> m (ReturnCode, Maybe (Double, Map v Double, [RowValue v c]))
-glpSolve' opts = get >>= liftIO . glpSolveAll opts
-
-{-# SPECIALIZE writeLPToFile :: (Ord v, Show v, Real c) => FilePath -> LPT v c IO () #-}
--- | Writes the current linear program to the specified file in CPLEX LP format. 
--- (This is a binding to GLPK, not a Haskell implementation of CPLEX.)
-writeLPToFile :: (Ord v, Show v, Real c, MonadState (LP v c) m, MonadIO m) =>
-	FilePath -> m ()
-writeLPToFile file = get >>= liftIO . writeLP file 
-
-{-# SPECIALIZE readLPFromFile :: (Ord v, Read v, Fractional c) => FilePath -> LPT v c IO () #-}
--- | Reads a linear program from the specified file in CPLEX LP format, overwriting
--- the current linear program.  Uses 'read' and 'realToFrac' to translate to the specified type.
--- Warning: this may not work on all files written using 'writeLPToFile', since variable names
--- may be changed.
--- (This is a binding to GLPK, not a Haskell implementation of CPLEX.)
-readLPFromFile :: (Ord v, Read v, Fractional c, MonadState (LP v c) m, MonadIO m) =>
-	FilePath -> m ()
-readLPFromFile = put <=< liftIO . readLP
-
-{-# SPECIALIZE readLPFromFile :: FilePath -> LPT String Double IO () #-}
--- | Reads a linear program from the specified file in CPLEX LP format, overwriting
--- the current linear program.  (This is a binding to GLPK, not a Haskell implementation of CPLEX.)
-readLPFromFile' :: (MonadState (LP String Double) m, MonadIO m) =>
-	FilePath -> m ()
-readLPFromFile' = put <=< liftIO . readLP'
diff --git a/Control/Monad/LPMonad/Internal.hs b/Control/Monad/LPMonad/Internal.hs
deleted file mode 100644
--- a/Control/Monad/LPMonad/Internal.hs
+++ /dev/null
@@ -1,245 +0,0 @@
-{-# LANGUAGE BangPatterns, FlexibleContexts, RecordWildCards #-}
-
-module Control.Monad.LPMonad.Internal (
---         module Data.LinearProgram.Common,
-        -- * Monad definitions
-        LPM,
-        LPT,
-        runLPM,
-        runLPT,
-        execLPM,
-        execLPT,
-        evalLPM,
-        evalLPT,
-        -- * Constructing the LP
-        -- ** Objective configuration
-        setDirection,
-        setObjective,
-        addObjective,
-        addWeightedObjective,
-        -- ** Two-function constraints
-        leq,
-        equal,
-        geq,
-        leq',
-        equal',
-        geq',
-        -- ** One-function constraints
-        leqTo,
-        equalTo,
-        geqTo,
-        constrain,
-        leqTo',
-        equalTo',
-        geqTo',
-        constrain',
-        -- ** Variable constraints
-        varLeq,
-        varEq,
-        varGeq,
-        varBds,
-        setVarBounds,
-        setVarKind,
---         newVariables,
---         newVariables'
-        ) where
-
-import Control.Monad.State.Strict
-import Control.Monad.Identity
-
-import Data.Map
-
-import Data.LinearProgram.Common
-
--- | A simple monad for constructing linear programs.  This library is intended to be able to link to 
--- a variety of different linear programming implementations.
-type LPM v c = LPT v c Identity
-
--- | A simple monad transformer for constructing linear programs in an arbitrary monad.
-type LPT v c = StateT (LP v c)
-
-runLPM :: (Ord v, Group c) => LPM v c a -> (a, LP v c)
-runLPM = runIdentity . runLPT
-
-runLPT :: (Ord v, Group c) => LPT v c m a -> m (a, LP v c)
-runLPT m = runStateT m (LP Max zero [] mempty mempty)
-
--- | Constructs a linear programming problem.
-execLPM :: (Ord v, Group c) => LPM v c a -> LP v c
-execLPM = runIdentity . execLPT
-
--- | Constructs a linear programming problem in the specified monad.
-execLPT :: (Ord v, Group c, Monad m) => LPT v c m a -> m (LP v c)
-execLPT = liftM snd . runLPT
-
--- | Runs the specified operation in the linear programming monad.
-evalLPM :: (Ord v, Group c) => LPM v c a -> a
-evalLPM = runIdentity . evalLPT
-
--- | Runs the specified operation in the linear programming monad transformer.
-evalLPT :: (Ord v, Group c, Monad m) => LPT v c m a -> m a
-evalLPT = liftM fst . runLPT
-
--- | Sets the optimization direction of the linear program: maximization or minimization.
-{-# SPECIALIZE setDirection :: Direction -> LPM v c (), Monad m => Direction -> LPT v c m () #-}
-setDirection :: (MonadState (LP v c) m) => Direction -> m ()
-setDirection dir = modify (\ lp -> lp{direction = dir})
-
-{-# SPECIALIZE equal :: (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 () #-}
-{-# SPECIALIZE leq :: (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 () #-}
-{-# 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.  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
-geq = flip leq
-
-{-# SPECIALIZE equal' :: (Ord v, Group c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
-        (Ord v, Group c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
-{-# SPECIALIZE geq' :: (Ord v, Group c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
-        (Ord v, Group c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
-{-# SPECIALIZE leq' :: (Ord v, Group c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
-        (Ord v, Group c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
--- | Specifies the relationship between two functions in the variables, with a label on the constraint.
-equal', leq', geq' :: (Ord v, Group c, MonadState (LP v c) m) => String -> LinFunc v c -> LinFunc v c -> m ()
-equal' lab f g = equalTo' lab (f ^-^ g) zero
-leq' lab f g = leqTo' lab (f ^-^ g) zero
-geq' = flip . leq'
-
-{-# SPECIALIZE equalTo :: LinFunc v c -> c -> LPM v c (), Monad m => LinFunc v c -> c -> LPT v c m () #-}
-{-# SPECIALIZE geqTo :: LinFunc v c -> c -> LPM v c (), Monad m => LinFunc v c -> c -> LPT v c m () #-}
-{-# SPECIALIZE leqTo :: LinFunc v c -> c -> LPM v c (), Monad m => LinFunc v c -> c -> LPT v c m () #-}
--- | Sets a constraint on a linear function in the variables.
-equalTo, leqTo, geqTo :: MonadState (LP v c) m => LinFunc v c -> c -> m ()
-equalTo f v = constrain f (Equ v)
-leqTo f v = constrain f (UBound v)
-geqTo f v = constrain f (LBound v)
-
-{-# SPECIALIZE equalTo' :: String -> LinFunc v c -> c -> LPM v c (),
-        Monad m => String -> LinFunc v c -> c -> LPT v c m () #-}
-{-# SPECIALIZE geqTo' :: String -> LinFunc v c -> c -> LPM v c (),
-        Monad m => String -> LinFunc v c -> c -> LPT v c m () #-}
-{-# SPECIALIZE leqTo' :: String -> LinFunc v c -> c -> LPM v c (),
-        Monad m => String -> LinFunc v c -> c -> LPT v c m () #-}
--- | Sets a labeled constraint on a linear function in the variables.
-equalTo', leqTo', geqTo' :: MonadState (LP v c) m => String -> LinFunc v c -> c -> m ()
-equalTo' lab f v = constrain' lab f (Equ v)
-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 = () <$ 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 () #-}
-{-# SPECIALIZE varLeq :: (Ord v, Ord c) => v -> c -> LPM v c (),
-        (Ord v, Ord c, Monad m) => v -> c -> LPT v c m () #-}
-{-# SPECIALIZE varGeq :: (Ord v, Ord c) => v -> c -> LPM v c (),
-        (Ord v, Ord c, Monad m) => v -> c -> LPT v c m () #-}
--- | Sets a constraint on the value of a variable.  If you constrain a variable more than once,
--- the constraints will be combined.  If the constraints are mutually contradictory,
--- an error will be generated.  This is more efficient than adding an equivalent function constraint.
-varEq, varLeq, varGeq :: (Ord v, Ord c, MonadState (LP v c) m) => v -> c -> m ()
-varEq v c = setVarBounds v (Equ c)
-varLeq v c = setVarBounds v (UBound c)
-varGeq v c = setVarBounds v (LBound c)
-
-{-# SPECIALIZE varBds :: (Ord v, Ord c) => v -> c -> c -> LPM v c (),
-        (Ord v, Ord c, Monad m) => v -> c -> c -> LPT v c m () #-}
--- | Bounds the value of a variable on both sides.  If you constrain a variable more than once,
--- the constraints will be combined.  If the constraints are mutually contradictory,
--- an error will be generated.  This is more efficient than adding an equivalent function constraint.
-varBds :: (Ord v, Ord c, MonadState (LP v c) m) => v -> c -> c -> m ()
-varBds v l u = setVarBounds v (Bound l u)
-
-{-# SPECIALIZE constrain :: LinFunc v c -> Bounds c -> LPM v c (),
-        Monad m => LinFunc v c -> Bounds c -> LPT v c m () #-}
--- | The most general form of an unlabeled constraint.
-constrain :: MonadState (LP v c) m => LinFunc v c -> Bounds c -> m ()
-constrain f bds = modify addConstr where
-        addConstr lp@LP{..}
-                = lp{constraints = Constr Nothing f bds:constraints}
-
-{-# SPECIALIZE constrain' :: String -> LinFunc v c -> Bounds c -> LPM v c (),
-        Monad m => String -> LinFunc v c -> Bounds c -> LPT v c m () #-}
--- | The most general form of a labeled constraint.
-constrain' :: MonadState (LP v c) m => String -> LinFunc v c -> Bounds c -> m ()
-constrain' lab f bds = modify addConstr where
-        addConstr lp@LP{..}
-                = lp{constraints = Constr (Just lab) f bds:constraints}
-
-{-# SPECIALIZE setObjective :: LinFunc v c -> LPM v c (),
-        Monad m => LinFunc v c -> LPT v c m () #-}
--- | Sets the objective function, overwriting the previous objective function.
-setObjective :: MonadState (LP v c) m => LinFunc v c -> m ()
-setObjective obj = modify setObj where
-        setObj lp = lp{objective = obj}
-
-{-# SPECIALIZE addObjective :: (Ord v, Group c) => LinFunc v c -> LPM v c (),
-        (Ord v, Group c, Monad m) => LinFunc v c -> LPT v c m () #-}
--- | Adds this function to the objective function.
-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
--- @'addObjective' (wt '*^' obj)@.
-addWeightedObjective :: (Ord v, Module r c, MonadState (LP v c) m) => r -> LinFunc v c -> m ()
-addWeightedObjective wt obj = addObjective (wt *^ obj)
-
-{-# SPECIALIZE setVarBounds :: (Ord v, Ord c) => v -> Bounds c -> LPM v c (),
-        (Ord v, Ord c, Monad m) => v -> Bounds c -> LPT v c m () #-}
--- | The most general way to set constraints on a variable.
--- If you constrain a variable more than once, the constraints will be combined.
--- If you combine mutually contradictory constraints, an error will be generated.
--- This is more efficient than creating an equivalent function constraint.
-setVarBounds :: (Ord v, Ord c, MonadState (LP v c) m) => v -> Bounds c -> m ()
-setVarBounds var bds = modify addBds where
-        addBds lp@LP{..} = lp{varBounds = insertWith mappend var bds varBounds}
-
-{-# SPECIALIZE setVarKind :: Ord v => v -> VarKind -> LPM v c (),
-        (Ord v, Monad m) => v -> VarKind -> LPT v c m () #-}
--- | Sets the kind ('type') of a variable.  See 'VarKind'.
-setVarKind :: (Ord v, MonadState (LP v c) m) => v -> VarKind -> m ()
-setVarKind v k = modify setK where
-        setK lp@LP{..} = lp{varTypes = insertWith mappend v k varTypes}
diff --git a/Control/Monad/LPMonad/Supply.hs b/Control/Monad/LPMonad/Supply.hs
deleted file mode 100644
--- a/Control/Monad/LPMonad/Supply.hs
+++ /dev/null
@@ -1,44 +0,0 @@
-{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}
-
-module Control.Monad.LPMonad.Supply (module Control.Monad.LPMonad.Supply.Class, Var(..), VSupply, VSupplyT, runVSupply, runVSupplyT) where
-
-import Control.Monad.Identity
-import Control.Monad.Trans
-import Control.Monad.State.Strict
-import Control.Monad.RWS.Class
-import Control.Monad.Cont.Class
-import Control.Monad.Error.Class
-import Control.Applicative
-import Control.Monad.LPMonad.Supply.Class
-
--- | A type suitable for use as a linear program variable.
-newtype Var = Var {varId :: Int} deriving (Eq, Ord, Enum)
-
--- | A monad capable of supplying unique variables.
-type VSupply = VSupplyT Identity
-
-runVSupply :: VSupply a -> a
-runVSupply = runIdentity . runVSupplyT
-
--- | A monad transformer capable of supplying unique variables.
-newtype VSupplyT m a = VSupplyT (StateT Var m a) deriving (Functor, Applicative, Monad, Alternative, MonadPlus, MonadTrans, MonadReader r, MonadWriter w, MonadCont,
-        MonadIO, MonadFix, MonadError e)
-
-runVSupplyT :: Monad m => VSupplyT m a -> m a
-runVSupplyT (VSupplyT m) = evalStateT m (Var 0)
-
-instance Show Var where
-        show (Var x) = "x_" ++ show x
-
-instance Read Var where
-        readsPrec _ ('x':'_':xs) = [(Var x, s') | (x, s') <- reads xs]
-        readsPrec _ _ = []
-
-instance MonadState s m => MonadState s (VSupplyT m) where
-        get = lift get
-        put = lift . put
-
-instance Monad m => MonadSupply Var (VSupplyT m) where
-        {-# SPECIALIZE instance MonadSupply Var VSupply #-}
-        supplyNew = VSupplyT $ StateT $ \ v -> return (v, succ v)
-        supplyN n = VSupplyT $ StateT $ \ (Var x) -> return (map Var [x..x+n-1], Var (x + n))
diff --git a/Control/Monad/LPMonad/Supply/Class.hs b/Control/Monad/LPMonad/Supply/Class.hs
deleted file mode 100644
--- a/Control/Monad/LPMonad/Supply/Class.hs
+++ /dev/null
@@ -1,52 +0,0 @@
-{-# LANGUAGE UndecidableInstances, FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies #-}
-module Control.Monad.LPMonad.Supply.Class where
-
-import Control.Monad
-
-import Control.Monad.State.Strict
-import Control.Monad.Reader
-import Control.Monad.Error
-import qualified Control.Monad.Writer.Lazy as WL
-import qualified Control.Monad.Writer.Strict as WS
-import qualified Control.Monad.State.Lazy as SL
-import Control.Monad.Cont
-
-import Data.Monoid
-
--- | A class implemented by monads that can supply values of type @s@.  Minimal implementation: 'supplyNew' or 'supplyN'.
-class Monad m => MonadSupply s m | m -> s where
-	-- | Supply a new value of type @s@.
-	supplyNew :: m s
-	-- | Supply @n@ values of type @s@.
-	supplyN :: Int -> m [s]
-	
-	supplyNew = liftM head (supplyN 1)
-	supplyN n = replicateM n supplyNew
-
-instance MonadSupply x m => MonadSupply x (StateT s m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
-
-instance MonadSupply x m => MonadSupply x (ReaderT r m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
-
-instance (Error e, MonadSupply x m) => MonadSupply x (ErrorT e m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
-
-instance (MonadSupply x m, Monoid w) => MonadSupply x (WL.WriterT w m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
-
-instance (MonadSupply x m, Monoid w) => MonadSupply x (WS.WriterT w m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
-
-instance MonadSupply x m => MonadSupply x (ContT r m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
-
-instance MonadSupply x m => MonadSupply x (SL.StateT s m) where
-	supplyNew = lift supplyNew
-	supplyN = lift . supplyN
diff --git a/Data/Algebra.hs b/Data/Algebra.hs
deleted file mode 100644
--- a/Data/Algebra.hs
+++ /dev/null
@@ -1,35 +0,0 @@
--- | Common library for algebraic structures.  Has the advantage of automatically inferring lots of useful structure, especially
--- in the writing of linear programs.  For example, here are several ways of writing @3 x - 4 y + z@:
--- 
--- > gsum [3 *& x, (-4) *^ var y, var z]
--- > linCombination [(3, x), (-4, y), (1, z)]
--- > 3 *& x ^-^ 4 *& y ^+^ var z
--- 
--- In addition, if we have two functions @f@ and @g@, we can construct linear combinations of those functions, using 
--- exactly the same syntax.  Moreover, we can multiply functions with 'Double' coefficients by 'Rational' values successfully.
--- This module is intended to offer as much generality as possible without getting in your way.
-module Data.Algebra (
-	-- * Algebraic structures
-	Group(..),
-	Ring(..),
-	Field(..),
-	Module(..),
-	VectorSpace(..),
-	Poly,
-	varPoly,
-	GroupRing,
-	LinFunc,
-	-- * Algebraic functions
-	gsum,
-	combination,
-	evalPoly,
-	-- ** Specialized methods on linear functions
-	var,
-	varSum,
-	(*&),
-	linCombination) where
-
-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
deleted file mode 100644
--- a/Data/Algebra/Field.hs
+++ /dev/null
@@ -1,24 +0,0 @@
-{-# 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/Group.hs b/Data/Algebra/Group.hs
deleted file mode 100644
--- a/Data/Algebra/Group.hs
+++ /dev/null
@@ -1,102 +0,0 @@
-{-# LANGUAGE TypeSynonymInstances #-}
-module Data.Algebra.Group where
-
-import Control.Applicative
-import qualified Data.Map as M
-import qualified Data.IntMap as IM
-import Data.Ratio
-
-type Poly = []
-
-infixr 4 ^+^
-infixr 4 ^-^
-
--- | The algebraic structure of a group.  Written additively.  Required functions: 'zero' and ('^-^' or ('^+^' and 'neg')).
-class Group g where
-	zero :: g
-	(^+^) :: g -> g -> g
-	(^-^) :: g -> g -> g
-	neg :: g -> g
-	
-	a ^+^ b = a ^-^ neg b
-	a ^-^ b = a ^+^ neg b
-	neg a = zero ^-^ a
-
-instance Group Bool where
-	zero = False
-	(^+^) = (/=)
-	(^-^) = (/=)
-	neg = id
-
-instance Group Int where
-	zero = 0
-	(^+^) = (+)
-	(^-^) = (-)
-	neg = negate
-
-instance Group Integer where
-	zero = 0
-	(^+^) = (+)
-	(^-^) = (-)
-	neg = negate
-
-instance Group Double where
-	zero = 0
-	(^+^) = (+)
-	(^-^) = (-)
-	neg = negate
-
-instance Integral a => Group (Ratio a) where
-	{-# SPECIALIZE instance Group Rational #-}
-	zero = 0
-	(^+^) = (+)
-	(^-^) = (-)
-	neg = negate
-
-instance Group g => Group (a -> g) where
-	zero = const zero
-	(^+^) = liftA2 (^+^)
-	(^-^) = liftA2 (^-^)
-	neg = fmap neg
-
-instance (Ord k, Group g) => Group (M.Map k g) where
-	zero = M.empty
-	(^+^) = M.unionWith (^+^)
-	neg = fmap neg
-
-instance Group g => Group (IM.IntMap g) where
-	zero = IM.empty
-	(^+^) = IM.unionWith (^+^)
-	neg = fmap neg
-
-instance Group g => Group (Poly g) where
-	zero = []
-	[] ^+^ p = p
-	p ^+^ [] = p
-	(a:as) ^+^ (b:bs) = (a ^+^ b):(as ^+^ bs)
-
-instance (Group g1, Group g2) => Group (g1, g2) where
-	{-# SPECIALIZE instance Group g => Group (g, g) #-}
-	zero = (zero, zero)
-	(x1, y1) ^+^ (x2, y2) = (x1 ^+^ x2, y1 ^+^ y2)
-	(x1, y1) ^-^ (x2, y2) = (x1 ^-^ x2, y1 ^-^ y2)
-	neg (x, y) = (neg x, neg y)
-
-instance (Group g1, Group g2, Group g3) => Group (g1, g2, g3) where
-	{-# SPECIALIZE instance Group g => Group (g, g, g) #-}
-	zero = (zero, zero, zero)
-	(x1, y1, z1) ^+^ (x2, y2, z2) = (x1 ^+^ x2, y1 ^+^ y2, z1 ^+^ z2)
-	(x1, y1, z1) ^-^ (x2, y2, z2) = (x1 ^-^ x2, y1 ^-^ y2, z1 ^-^ z2)
-	neg (x, y, z) = (neg x, neg y, neg z)
-
-instance (Group g1, Group g2, Group g3, Group g4) => Group (g1, g2, g3, g4) where
-	{-# SPECIALIZE instance Group g => Group (g, g, g, g) #-}
-	zero = (zero, zero, zero, zero)
-	(x1, y1, z1, w1) ^+^ (x2, y2, z2, w2) = (x1 ^+^ x2, y1 ^+^ y2, z1 ^+^ z2, w1 ^+^ w2)
-	(x1, y1, z1, w1) ^-^ (x2, y2, z2, w2) = (x1 ^-^ x2, y1 ^-^ y2, z1 ^-^ z2, w1 ^-^ w2)
-	neg (x, y, z, w) = (neg x, neg y, neg z, neg w)
-
-{-# INLINE gsum #-}
--- | Does a summation over the elements of a group.
-gsum :: Group g => [g] -> g
-gsum = foldr (^+^) zero
diff --git a/Data/Algebra/Module.hs b/Data/Algebra/Module.hs
deleted file mode 100644
--- a/Data/Algebra/Module.hs
+++ /dev/null
@@ -1,109 +0,0 @@
-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, IncoherentInstances, TypeSynonymInstances #-}
-
-module Data.Algebra.Module where
-
-import Data.Ratio
-import qualified Data.Map as M
-import qualified Data.IntMap as IM
-
-import Data.Algebra.Group
-import Data.Algebra.Ring
-
--- | The algebraic structure of a module.  A vector space is a module with coefficients in a field.
-class (Ring r, Group m) => Module r m where
-	(*^) :: r -> m -> m
-
-instance Module Int Int where
-	(*^) = (*)
-
-instance Module Integer Integer where
-	(*^) = (*)
-
-instance Module Int Integer where
-	(*^) = (*) . fromIntegral
-
-instance Integral a => Module Int (Ratio a) where
-	{-# SPECIALIZE instance Module Int Rational #-}
-	(*^) = (*) . fromIntegral
-
-instance Integral a => Module Integer (Ratio a) where
-	{-# SPECIALIZE instance Module Integer Rational #-}
-	(*^) = (*) . fromIntegral
-
-instance Integral a => Module (Ratio a) (Ratio a) where
-	{-# SPECIALIZE instance Module Rational Rational #-}
-	(*^) = (*)
-
-instance Module Int Double where
-	(*^) = (*) . fromIntegral
-
-instance Module Integer Double where
-	(*^) = (*) . fromIntegral
-
-instance Integral a => Module (Ratio a) Double where
-	{-# SPECIALIZE instance Module Rational Double #-}
-	(*^) = (*) . realToFrac
-
-instance Module Double Double where
-	(*^) = (*)
-
-instance (Ord g, Group g, Ring r) => Module (GroupRing r g) (GroupRing r g) where
-	(*^) = (*#)
-
-instance Module r m => Module r (a -> m) where
-	(*^) = fmap . (*^)
-
-instance (Ord k, Module r m) => Module r (M.Map k m) where
-	(*^) = fmap . (*^)
-
-instance Module r m => Module r (IM.IntMap m) where
-	(*^) = fmap . (*^)
-
-instance (Module r m1, Module r m2) => Module r (m1, m2) where
-	{-# SPECIALIZE instance Module r m => Module r (m, m) #-}
-	r *^ (a, b) = (r *^ a, r *^ b)
-
-instance (Module r m1, Module r m2, Module r m3) => Module r (m1, m2, m3) where
-	{-# SPECIALIZE instance Module r m => Module r (m, m, m) #-}
-	r *^ (a, b, c) = (r *^ a, r *^ b, r *^ c)
-
-instance (Module r m1, Module r m2, Module r m3, Module r m4) => Module r (m1, m2, m3, m4) where
-	{-# SPECIALIZE instance Module r m => Module r (m, m, m, m) #-}
-	r *^ (a, b, c, d) = (r *^ a, r *^ b, r *^ c, r *^ d)
-
--- | @'LinFunc' v c@ is a linear combination of variables of type @v@ with coefficients
--- from @c@.  Formally, this is the free @c@-module on @v@.  
-type LinFunc = M.Map
-
--- | Given a variable @v@, returns the function equivalent to @v@.
-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, Ring c) => c -> v -> LinFunc v c
-c *& v = M.singleton v c
-
--- | Equivalent to @'vsum' . 'map' 'var'@.
-varSum :: (Ord v, Ring c) => [v] -> LinFunc v c
-varSum vs = M.fromList [(v, one) | v <- vs]
-
--- | Given a collection of vectors and scaling coefficients, returns this
--- linear combination.
-combination :: Module r m => [(r, m)] -> m
-combination xs = gsum [r *^ m | (r, m) <- xs]
-
-{-# INLINE linCombination #-}
--- | Given a set of basic variables and coefficients, returns the linear combination obtained
--- by summing.
-linCombination :: (Ord v, Num r) => [(r, v)] -> LinFunc v r
-linCombination xs = M.fromListWith (+) [(v, r) | (r, v) <- xs]
-
--- | Substitution into a polynomial.
-evalPoly :: (Module r m, Ring m) => Poly r -> m -> m
-evalPoly f x = foldr (\ c z -> (c *^ one) ^+^ (x *# z)) zero f
-
-{-# RULES
-	"zero/*^" forall m . zero *^ m = zero;
-	"*^/zero" forall r . r *^ zero = zero;
-	"one/*^" forall m . one *^ m = m;
-	#-}
diff --git a/Data/Algebra/Ring.hs b/Data/Algebra/Ring.hs
deleted file mode 100644
--- a/Data/Algebra/Ring.hs
+++ /dev/null
@@ -1,61 +0,0 @@
-{-# LANGUAGE TypeSynonymInstances #-}
-module Data.Algebra.Ring where
-
-import Control.Applicative
-
-import Data.Ratio
-import qualified Data.Map as M
-
-import Data.Algebra.Group
-
-infixr 6 *#
-
--- | A way of forming a ring from functions.  See <http://en.wikipedia.org/wiki/Group_ring>.
-type GroupRing r g = M.Map g r
-
--- | The algebraic structure of a unital ring.  Assumes that the additive operation forms an abelian group,
--- that the multiplication operation forms a group, and that multiplication distributes.
-class Group r => Ring r where
-	one :: r
-	(*#) :: r -> r -> r
-
-instance Ring Bool where
-	one = True
-	(*#) = (&&)
-
-instance Ring Int where
-	one = 1
-	(*#) = (*)
-
-instance Ring Integer where
-	one = 1
-	(*#) = (*)
-
-instance Ring Double where
-	one = 1
-	(*#) = (*)
-
-instance Integral a => Ring (Ratio a) where
-	{-# SPECIALIZE instance Ring Rational #-}
-	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]
-
--- | Returns the polynomial @p(x) = x@.
-varPoly :: Ring r => Poly r
-varPoly = [zero, one]
diff --git a/Data/LinearProgram.hs b/Data/LinearProgram.hs
deleted file mode 100644
--- a/Data/LinearProgram.hs
+++ /dev/null
@@ -1,8 +0,0 @@
-module Data.LinearProgram (
-	module Data.LinearProgram.Common,
-	module Data.LinearProgram.GLPK,
-	module Control.Monad.LPMonad) where
-
-import Data.LinearProgram.GLPK
-import Data.LinearProgram.Common
-import Control.Monad.LPMonad
diff --git a/Data/LinearProgram/Common.hs b/Data/LinearProgram/Common.hs
deleted file mode 100644
--- a/Data/LinearProgram/Common.hs
+++ /dev/null
@@ -1,19 +0,0 @@
--- | Contains sufficient tools to represent linear programming problems in Haskell.  In the future, if linkings to other
--- linear programming libraries are made, this will be common to them all.
-module Data.LinearProgram.Common (
-	module Data.LinearProgram.Spec,
-	module Data.Algebra,
-	module Data.LinearProgram.Types) where
-
-import Data.LinearProgram.Spec
-import Data.Algebra
-import Data.LinearProgram.Types
-
-import Data.Map
-import GHC.Exts (build)
-
-{-# RULES
-	"assocs" assocs = \ m -> build (\ c n -> foldWithKey (curry c) n m);
-	"elems" elems = \ m -> build (\ c n -> foldWithKey (const c) n m);
-	"keys" keys = \ m -> build (\ c n -> foldWithKey (\ k _ -> c k) n m);
-	#-}
diff --git a/Data/LinearProgram/GLPK.hs b/Data/LinearProgram/GLPK.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK.hs
+++ /dev/null
@@ -1,8 +0,0 @@
-module Data.LinearProgram.GLPK (
--- 	module Data.LinearProgram.LPMonad,
-	module Data.LinearProgram.GLPK.Solver,
-	module Data.LinearProgram.GLPK.IO) where
-	
--- import Data.LinearProgram.LPMonad
-import Data.LinearProgram.GLPK.Solver
-import Data.LinearProgram.GLPK.IO
diff --git a/Data/LinearProgram/GLPK/Common.hs b/Data/LinearProgram/GLPK/Common.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK/Common.hs
+++ /dev/null
@@ -1,13 +0,0 @@
-module Data.LinearProgram.GLPK.Common (
-	module Data.LinearProgram.GLPK.Internal,
-	module Data.LinearProgram.GLPK.Types,
-	module Foreign.Ptr,
-	module Foreign.C,
-	module Foreign.Marshal.Array) where
-
-import Data.LinearProgram.GLPK.Internal
-import Data.LinearProgram.GLPK.Types
-
-import Foreign.Ptr
-import Foreign.C
-import Foreign.Marshal.Array
diff --git a/Data/LinearProgram/GLPK/IO.hs b/Data/LinearProgram/GLPK/IO.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK/IO.hs
+++ /dev/null
@@ -1,21 +0,0 @@
--- | Bindings to the file I/O functions from GLPK, on the CPLEX LP file format.
-module Data.LinearProgram.GLPK.IO where
-
-import Data.LinearProgram.Common
-
-import Data.LinearProgram.GLPK.Common
-import Data.LinearProgram.GLPK.IO.Internal
-
-{-# SPECIALIZE readLP :: (Ord v, Read v) => FilePath -> IO (LP v Double) #-}
--- | Read a linear program from a file in CPLEX LP format.  Warning: this will not necessarily succeed
--- on all files generated by 'writeLP', as variable names may be changed.
-readLP :: (Ord v, Read v, Fractional c) => FilePath -> IO (LP v c)
-readLP = fmap (mapVals realToFrac . mapVars read) . readLP'
-
--- | Read a linear program from a file in CPLEX LP format.
-readLP' :: FilePath -> IO (LP String Double)
-readLP' = runGLPK . readGLPLP
-
--- | Write a linear program to a file in CPLEX LP format.
-writeLP :: (Ord v, Show v, Real c) => FilePath -> LP v c -> IO ()
-writeLP file = runGLPK . writeGLPLP file
diff --git a/Data/LinearProgram/GLPK/IO/Internal.hs b/Data/LinearProgram/GLPK/IO/Internal.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK/IO/Internal.hs
+++ /dev/null
@@ -1,134 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
-module Data.LinearProgram.GLPK.IO.Internal (readGLPLP, writeGLPLP) where
-
-import Control.Monad
-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
-import Control.Monad.LPMonad.Internal
-
-foreign import ccall unsafe "c_glp_write_lp" glpWriteLP :: Ptr GlpProb -> CString -> IO ()
-foreign import ccall unsafe "c_glp_read_lp" glpReadLP :: Ptr GlpProb -> CString -> IO ()
-foreign import ccall unsafe "c_glp_set_col_name" glpSetColName :: Ptr GlpProb -> CInt -> CString -> IO ()
-foreign import ccall unsafe "c_glp_set_row_name" glpSetRowName :: Ptr GlpProb -> CInt -> CString -> IO ()
-foreign import ccall unsafe "c_glp_get_obj_dir" glpGetObjDir :: Ptr GlpProb -> IO CInt
-foreign import ccall unsafe "c_glp_get_num_rows" glpGetNumRows :: Ptr GlpProb -> IO CInt
-foreign import ccall unsafe "c_glp_get_num_cols" glpGetNumCols :: Ptr GlpProb -> IO CInt
-foreign import ccall unsafe "c_glp_get_row_name" glpGetRowName :: Ptr GlpProb -> CInt -> IO CString
-foreign import ccall unsafe "c_glp_get_col_name" glpGetColName :: Ptr GlpProb -> CInt -> IO CString
-foreign import ccall unsafe "c_glp_get_col_kind" glpGetColKind :: Ptr GlpProb -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_get_row_type" glpGetRowType :: Ptr GlpProb -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_get_col_type" glpGetColType :: Ptr GlpProb -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_get_row_lb" glpGetRowLb :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_get_col_lb" glpGetColLb :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_get_row_ub" glpGetRowUb :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_get_col_ub" glpGetColUb :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_get_obj_coef" glpGetObjCoef :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_get_mat_row" glpGetMatRow :: Ptr GlpProb -> CInt -> Ptr CInt -> Ptr CDouble -> IO CInt
-
-writeLP :: FilePath -> GLPK ()
-writeLP file = GLP $ withCString file . glpWriteLP
-
-readLP :: FilePath -> GLPK ()
-readLP file = GLP $ withCString file . glpReadLP
-
-getDir :: GLPK Direction
-getDir = liftM (toEnum . subtract 1 . fromIntegral) $ GLP glpGetObjDir
-
-getRowName, getColName :: Int -> GLPK (Maybe String)
-getRowName i = GLP $ peekCAString' <=< flip glpGetRowName (fromIntegral i)
-getColName i = GLP $ peekCAString' <=< flip glpGetColName (fromIntegral i)
-
-peekCAString' :: CString -> IO (Maybe String)
-peekCAString' str
-	| str == nullPtr	= return Nothing
-	| otherwise		= liftM Just $ peekCAString str
-
-getNumRows, getNumCols :: GLPK Int
-getNumRows = liftM fromIntegral $ GLP glpGetNumRows
-getNumCols = liftM fromIntegral $ GLP glpGetNumCols
-
-rowBounds, colBounds :: Int -> GLPK (Bounds Double)
-rowBounds = loadBounds (getCDouble glpGetRowLb) (getCDouble glpGetRowUb) (getCInt glpGetRowType)
-colBounds = loadBounds (getCDouble glpGetColLb) (getCDouble glpGetColUb) (getCInt glpGetColType)
-
-colKind :: Int -> GLPK VarKind
-colKind = liftM (toEnum . subtract 1) . getCInt glpGetColKind
-
-getCInt :: (Ptr GlpProb -> CInt -> IO CInt) -> Int -> GLPK Int
-getCInt f i = GLP $ \ lp -> liftM fromIntegral $ f lp (fromIntegral i)
-
-getCDouble :: (Ptr GlpProb -> CInt -> IO CDouble) -> Int -> GLPK Double
-getCDouble f i = GLP $ \ lp -> liftM realToFrac $ f lp (fromIntegral i)
-
-setRowName :: Int -> String -> GLPK ()
-setRowName i nam = GLP $ withCString nam . flip glpSetRowName (fromIntegral i)
-
-setColName :: Int -> String -> GLPK ()
-setColName i nam = GLP $ withCString nam . flip glpSetColName (fromIntegral i)
-
-loadBounds :: (Int -> GLPK Double) -> (Int -> GLPK Double) ->
-	(Int -> GLPK Int) -> Int -> GLPK (Bounds Double)
-loadBounds lb ub tp i = do
-	typ <- tp i
-	case typ of
-		1	-> return Free
-		2	-> liftM LBound (lb i)
-		3	-> liftM UBound (ub i)
-		4	-> liftM2 Bound (lb i) (ub i)
-		_	-> liftM Equ (lb i)
-		
-getObjCoef :: Int -> GLPK Double
-getObjCoef = getCDouble glpGetObjCoef
-
-getRows :: GLPK [(Int, [(Int, Double)])]
-getRows = do	n <- getNumRows
-		m <- getNumCols
-		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 $ mapM (peekElemOff ixs) [1..k]
-			coefsL <- liftIO $ mapM (peekElemOff ixs) [1..k]
-			return (i, zip (map fromIntegral ixsL) (map realToFrac coefsL))
-			| i <- [1..n]]
-
-readGLPLP :: FilePath -> GLPK (LP String Double)
-readGLPLP file = execLPT $ do
-	lift $ readLP file
-	setDirection =<< lift getDir
-	nCols <- lift getNumCols
-	names <- lift $ liftM fromList $ mapM (\ i -> do
-		Just name <- getColName i
-		return (i, name)) [1..nCols]
-	sequence_ [do
-		bds <- lift $ colBounds i
-		kind <- lift $ colKind i
-		setVarBounds name bds
-		setVarKind name kind
-		return (i, name)
-			| (i, name) <- assocs names]
-	rowContents <- lift getRows
-	sequence_ [do
-		bds <- lift $ rowBounds i
-		name <- lift $ getRowName i
-		maybe constrain constrain' name 
-			(linCombination [(v, names ! j) | (j, v) <- row]) bds
-			| (i, row) <- rowContents]
-	obj <- lift $ sequence [do
-		c <- getObjCoef i
-		return (name, c) | (i, name) <- assocs names]
-	setObjective (fromList (filter ((/= 0) . snd) obj))
-
-writeGLPLP :: (Show v, Ord v, Real c) => FilePath -> LP v c -> GLPK ()
-writeGLPLP file lp = do
-	vars <- writeProblem lp
-	sequence_ [setColName i (show v) | (v, i) <- assocs vars]
-	sequence_ [setRowName i lab | (i, Constr (Just lab) _ _) <- zip [1..] (constraints lp)]
-	writeLP file
diff --git a/Data/LinearProgram/GLPK/Internal.hs b/Data/LinearProgram/GLPK/Internal.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK/Internal.hs
+++ /dev/null
@@ -1,173 +0,0 @@
-{-# LANGUAGE RecordWildCards, ScopedTypeVariables, ForeignFunctionInterface, BangPatterns #-}
-module Data.LinearProgram.GLPK.Internal (writeProblem, solveSimplex, mipSolve,
-	getObjVal, getRowPrim, getColPrim, mipObjVal, mipRowVal, mipColVal, getBadRay) where
-{-(writeProblem, addCols,
-	addRows, createIndex, findCol, findRow, getColPrim, getRowPrim, getObjVal,
-	mipColVal, mipRowVal, mipObjVal, mipSolve, setColBounds, setColKind, setColName, setMatRow,
-	setObjCoef, setObjectiveDirection, setRowBounds, setRowName, solveSimplex) where-}
-
-import Control.Monad
-
-import Foreign.Ptr
-import Foreign.C
-import Foreign.Marshal.Array
-
-import Data.Bits
-import Data.Map hiding (map)
--- import Data.Bounds
-import Data.LinearProgram.Common
-import Data.LinearProgram.GLPK.Types
-
--- foreign import ccall "c_glp_set_obj_name" glpSetObjName :: Ptr GlpProb -> CString -> IO ()
--- foreign import ccall unsafe "c_glp_set_obj_dir" glpSetObjDir :: Ptr GlpProb -> CInt -> IO ()
-foreign import ccall unsafe "c_glp_minimize" glpMinimize :: Ptr GlpProb -> IO ()
-foreign import ccall unsafe "c_glp_maximize" glpMaximize :: Ptr GlpProb -> IO ()
-foreign import ccall unsafe "c_glp_add_rows" glpAddRows :: Ptr GlpProb -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_add_cols" glpAddCols :: Ptr GlpProb -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_set_row_bnds" glpSetRowBnds :: Ptr GlpProb -> CInt -> CInt -> CDouble -> CDouble -> IO ()
-foreign import ccall unsafe "c_glp_set_col_bnds" glpSetColBnds :: Ptr GlpProb -> CInt -> CInt -> CDouble -> CDouble -> IO ()
-foreign import ccall unsafe "c_glp_set_obj_coef" glpSetObjCoef :: Ptr GlpProb -> CInt -> CDouble -> IO ()
-foreign import ccall unsafe "c_glp_set_mat_row" glpSetMatRow :: Ptr GlpProb -> CInt -> CInt -> Ptr CInt -> Ptr CDouble -> IO ()
--- foreign import ccall unsafe "c_glp_create_index" glpCreateIndex :: Ptr GlpProb -> IO ()
--- foreign import ccall unsafe "c_glp_find_row" glpFindRow :: Ptr GlpProb -> CString -> IO CInt
--- foreign import ccall unsafe "c_glp_find_col" glpFindCol :: Ptr GlpProb -> CString -> IO CInt
-foreign import ccall unsafe "c_glp_solve_simplex" glpSolveSimplex :: Ptr GlpProb -> CInt -> CInt -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_get_obj_val" glpGetObjVal :: Ptr GlpProb -> IO CDouble
-foreign import ccall unsafe "c_glp_get_row_prim" glpGetRowPrim :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_get_col_prim" glpGetColPrim :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_set_col_kind" glpSetColKind :: Ptr GlpProb -> CInt -> CInt -> IO ()
-foreign import ccall unsafe "c_glp_mip_solve" glpMipSolve :: 
-	Ptr GlpProb -> CInt -> CInt -> CInt -> CInt -> CInt -> CInt -> CInt -> CDouble -> CInt -> IO CInt
-foreign import ccall unsafe "c_glp_mip_obj_val" glpMIPObjVal :: Ptr GlpProb -> IO CDouble
-foreign import ccall unsafe "c_glp_mip_row_val" glpMIPRowVal :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_mip_col_val" glpMIPColVal :: Ptr GlpProb -> CInt -> IO CDouble
-foreign import ccall unsafe "c_glp_set_row_name" glpSetRowName :: Ptr GlpProb -> CInt -> CString -> IO ()
-foreign import ccall unsafe "c_glp_get_bad_ray" glpGetBadRay :: Ptr GlpProb -> IO CInt
-
-setObjectiveDirection :: Direction -> GLPK ()
-setObjectiveDirection dir = GLP $ case dir of
-	Min	-> glpMinimize
-	Max	-> glpMaximize
-
-getBadRay :: GLPK (Maybe Int)
-getBadRay = liftM (\ x -> guard (x /= 0) >> return (fromIntegral x)) $ GLP glpGetBadRay
-
-addRows :: Int -> GLPK Int
-addRows n = GLP $ liftM fromIntegral . flip glpAddRows (fromIntegral n)
-
-addCols :: Int -> GLPK Int
-addCols n = GLP $ liftM fromIntegral . flip glpAddCols (fromIntegral n)
-
-setRowBounds :: Real a => Int -> Bounds a -> GLPK ()
-setRowBounds i bds = GLP $ \ lp -> onBounds (glpSetRowBnds lp (fromIntegral i)) bds
-
-setColBounds :: Real a => Int -> Bounds a -> GLPK ()
-setColBounds i bds = GLP $ \ lp -> onBounds (glpSetColBnds lp (fromIntegral i)) bds
-
-onBounds :: Real a => (CInt -> CDouble -> CDouble -> x) -> Bounds a -> x
-onBounds f bds = case bds of
-	Free		-> f 1 0 0
-	LBound a	-> f 2 (realToFrac a) 0
-	UBound a	-> f 3 0 (realToFrac a)
-	Bound a b	-> f 4 (realToFrac a) (realToFrac b)
-	Equ a		-> f 5 (realToFrac a) 0
-
-{-# SPECIALIZE setObjCoef :: Int -> Double -> GLPK (), Int -> Int -> GLPK () #-}
-setObjCoef :: Real a => Int -> a -> GLPK ()
-setObjCoef i v = GLP $ \ lp -> glpSetObjCoef lp (fromIntegral i) (realToFrac v)
-
-{-# SPECIALIZE setMatRow :: Int -> [(Int, Double)] -> GLPK (), Int -> [(Int, Int)] -> GLPK () #-}
-setMatRow :: Real a => Int -> [(Int, a)] -> GLPK ()
-setMatRow i row = GLP $ \ lp -> 
-	allocaArray (len+1) $ \ (ixs :: Ptr CInt) -> allocaArray (len+1) $ \ (coeffs :: Ptr CDouble) -> do
-		pokeArray ixs (0:map (fromIntegral . fst) row)
-		pokeArray coeffs (0:map (realToFrac . snd) row)
-		glpSetMatRow lp (fromIntegral i) (fromIntegral len) ixs coeffs
-	where	len = length row
-
--- createIndex :: GLPK ()
--- createIndex = GLP glpCreateIndex
-
--- findRow :: String -> GLPK Int
--- findRow nam = GLP $ liftM fromIntegral . withCString nam . glpFindRow
-
--- findCol :: String -> GLPK Int
--- findCol nam = GLP $ liftM fromIntegral . withCString nam . glpFindCol
-
-solveSimplex :: MsgLev -> Int -> Bool -> GLPK ReturnCode
-solveSimplex msglev tmLim presolve = GLP $ \ lp -> liftM (toEnum . fromIntegral) $ glpSolveSimplex lp
-	(getMsgLev msglev)
-	tmLim'
-	(if presolve then 1 else 0)
-	where	tmLim' = fromIntegral (tmLim * 1000)
-
-getMsgLev :: MsgLev -> CInt
-getMsgLev = fromIntegral . fromEnum
-
-getObjVal :: GLPK Double
-getObjVal = liftM realToFrac $ GLP glpGetObjVal
-
-getRowPrim :: Int -> GLPK Double
-getRowPrim i = liftM realToFrac $ GLP (`glpGetRowPrim` fromIntegral i)
-
-getColPrim :: Int -> GLPK Double
-getColPrim i = liftM realToFrac $ GLP (`glpGetColPrim` fromIntegral i)
-
-setColKind :: Int -> VarKind -> GLPK ()
-setColKind i kind = GLP $ \ lp -> glpSetColKind lp (fromIntegral i) (fromIntegral $ 1 + fromEnum kind)
-
-mipSolve :: MsgLev -> BranchingTechnique -> BacktrackTechnique -> Preprocessing -> Bool ->
-	[Cuts] -> Double -> Int -> Bool -> GLPK ReturnCode
-mipSolve msglev brt btt pp fp cuts mipgap tmlim presol =
-		liftM (toEnum . fromIntegral) $ GLP $ \ lp -> glpMipSolve lp msglev'
-						brt' btt' pp' fp' tmlim' cuts' mipgap' presol'
-	where	!msglev' = getMsgLev msglev
-		!brt' = 1 + fromIntegral (fromEnum brt)
-		!btt' = 1 + fromIntegral (fromEnum btt)
-		!pp' = fromIntegral (fromEnum pp)
-		!fp' = fromIntegral (fromEnum fp)
-		!cuts' = (if GMI `elem` cuts then 1 else 0) .|.
-			(if MIR `elem` cuts then 2 else 0) .|.
-			(if Cov `elem` cuts then 4 else 0) .|.
-			(if Clq `elem` cuts then 8 else 0)
-		!mipgap' = realToFrac mipgap
-		!tmlim' = fromIntegral (1000 * tmlim)
-		!presol' = fromIntegral (fromEnum presol)
-
-mipObjVal :: GLPK Double
-mipObjVal = liftM realToFrac $ GLP glpMIPObjVal
-
-mipRowVal :: Int -> GLPK Double
-mipRowVal i = liftM realToFrac $ GLP (`glpMIPRowVal` fromIntegral i)
-
-mipColVal :: Int -> GLPK Double
-mipColVal i = liftM realToFrac $ GLP (`glpMIPColVal` fromIntegral i)
-
-setRowName :: Int -> String -> GLPK ()
-setRowName i nam = GLP $ withCString nam . flip glpSetRowName (fromIntegral i)
-
-{-# SPECIALIZE writeProblem :: Ord v => LP v Double -> GLPK (Map v Int),
-	Ord v => LP v Int -> GLPK (Map v Int) #-}
-writeProblem :: (Ord v, Real c) => LP v c -> GLPK (Map v Int)
-writeProblem LP{..} = do
-	setObjectiveDirection direction
-	i0 <- addCols nVars
-	let allVars' = fmap (i0 +) allVars
-	sequence_ [setObjCoef i v | (i, v) <- elems $ intersectionWith (,) allVars' objective]
-	j0 <- addRows (length constraints)
-	sequence_ [do	maybe (return ()) (setRowName j) lab
-			setMatRow j
-				[(i, v) | (i, v) <- elems (intersectionWith (,) allVars' f)]
-			setRowBounds j bnds
-				| (j, Constr lab f bnds) <- zip [j0..] constraints]
--- 	createIndex
-	sequence_ [setColBounds i bnds |
-			(i, bnds) <- elems $ intersectionWith (,) allVars' varBounds]
-	sequence_ [setColBounds i Free | i <- elems $ difference allVars' varBounds]
-	sequence_ [setColKind i knd |
-			(i, knd) <- elems $ intersectionWith (,) allVars' varTypes]
-	return allVars'
-	where	allVars0 = fmap (const ()) objective `union`
-			unions [fmap (const ()) f | Constr _ f _ <- constraints] `union`
-			fmap (const ()) varBounds `union` fmap (const ()) varTypes
-		(nVars, allVars) = mapAccum (\ n _ -> (n+1, n)) (0 :: Int) allVars0
diff --git a/Data/LinearProgram/GLPK/Solver.hs b/Data/LinearProgram/GLPK/Solver.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK/Solver.hs
+++ /dev/null
@@ -1,119 +0,0 @@
-{-# OPTIONS -funbox-strict-fields #-}
-{-# LANGUAGE TupleSections, RecordWildCards #-}
-
--- | Interface between the Haskell representation of a linear programming problem, a value of type 'LP', and
--- the GLPK solver.  The options available to the solver correspond naturally with GLPK's available options,
--- so to find the meaning of any particular option, consult the GLPK documentation.
--- 
--- The option of which solver to use -- the general LP solver, which solves a problem over the reals, or the 
--- MIP solver, which allows variables to be restricted to integers -- can be made by choosing the appropriate
--- constructor for 'GLPOpts'.
--- 
--- The marshalling from Haskell to C is specialized for 'Int's and 'Double's, so using those types in your
--- linear program is recommended.
-module Data.LinearProgram.GLPK.Solver (
-	-- * Solver options
-	GLPOpts(..),
-	simplexDefaults, 
-	mipDefaults, 
-	-- * Running the solver
-	glpSolveVars,
-	RowValue(..),
-	glpSolveAll,
-	-- * GLPK enumerations
-	ReturnCode(..),
-	MsgLev(..), 
-	BranchingTechnique(..),
-	BacktrackTechnique(..), 
-	Preprocessing(..), 
-	Cuts(..)) where 
-
-import Control.Monad
-
-import Data.Map
-import Data.LinearProgram.Spec
-import Data.LinearProgram.GLPK.Common
-
--- | Options available for customizing GLPK operations.  This also determines
--- which kind of solving is performed -- relaxed LP, or MIP.
-data GLPOpts = SimplexOpts {msgLev :: MsgLev, tmLim :: !Int, presolve :: Bool} |
-	MipOpts {msgLev :: MsgLev, tmLim :: !Int, presolve :: Bool,
-		brTech :: BranchingTechnique, btTech :: BacktrackTechnique,
-		ppTech :: Preprocessing,
-		fpHeur :: Bool,
-		cuts :: [Cuts],
-		mipGap :: !Double}
-
-data RowValue v c = RowVal {row :: !(Constraint v c), rowVal :: !Double}
-
-simplexDefaults, mipDefaults :: GLPOpts
-simplexDefaults = SimplexOpts MsgOn 10000 True
-mipDefaults = MipOpts MsgOn 10000 True DrTom LocBound AllPre False [] 0.0
-
-{-# SPECIALIZE glpSolveVars :: Ord v => GLPOpts -> LP v Double -> IO (ReturnCode, Maybe (Double, Map v Double)),
-	Ord v => GLPOpts -> LP v Int -> IO (ReturnCode, Maybe (Double, Map v Double)) #-}
--- | Solves the linear or mixed integer programming problem.  Returns
--- the value of the objective function, and the values of the variables.
-glpSolveVars :: (Ord v, Real c) => GLPOpts -> LP v c -> IO (ReturnCode, Maybe (Double, Map v Double))
-glpSolveVars opts@SimplexOpts{} lp = runGLPK $ do
-	(code, vars) <- doGLP opts lp
-	liftM (code, ) $ maybe (return Nothing) ( \ vars -> do
-		obj <- getObjVal
-		vals <- sequence [do
-			val <- getColPrim i
-			return (v, val)
-				| (v, i) <- assocs vars]
-		return (Just (obj, fromDistinctAscList vals))) vars
-glpSolveVars opts@MipOpts{} lp = runGLPK $ do
-	(code, vars) <- doGLP opts lp
-	liftM (code, ) $ maybe (return Nothing) (\ vars -> do
-		obj <- mipObjVal
-		vals <- sequence [do
-			val <- mipColVal i
-			return (v, val)
-				| (v, i) <- assocs vars]
-		return (Just (obj, fromDistinctAscList vals))) vars
-
-{-# SPECIALIZE glpSolveAll :: 
-	Ord v => GLPOpts -> LP v Double -> IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v Double])),
-	Ord v => GLPOpts -> LP v Int -> IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v Int])) #-}
--- | Solves the linear or mixed integer programming problem.  Returns
--- the value of the objective function, the values of the variables,
--- and the values of any labeled rows.
-glpSolveAll :: (Ord v, Real c) => GLPOpts -> LP v c -> IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c]))
-glpSolveAll opts@SimplexOpts{} lp@LP{..} = runGLPK $ do
-	(code, vars) <- doGLP opts lp
-	liftM (code, ) $ maybe (return Nothing) (\ vars -> do
-		obj <- getObjVal
-		vals <- sequence [do
-			val <- getColPrim i
-			return (v, val)
-				| (v, i) <- assocs vars]
-		rows <- sequence [liftM (RowVal c) (getRowPrim i)
-					| (i, c) <- zip [1..] constraints]
-		return (Just (obj, fromDistinctAscList vals, rows))) vars
-glpSolveAll opts@MipOpts{} lp@LP{..} = runGLPK $ do
-	(code, vars) <- doGLP opts lp
-	liftM (code, ) $ maybe (return Nothing) (\ vars -> do
-		obj <- mipObjVal
-		vals <- sequence [do
-			val <- mipColVal i
-			return (v, val)
-				| (v, i) <- assocs vars]
-		rows <- sequence [liftM (RowVal c) (mipRowVal i)
-					| (i, c) <- zip [1..] constraints]
-		return (Just (obj, fromDistinctAscList vals, rows))) vars
-
-{-# SPECIALIZE doGLP :: Ord v => GLPOpts -> LP v Double -> GLPK (ReturnCode, Maybe (Map v Int)),
-	Ord v => GLPOpts -> LP v Int -> GLPK (ReturnCode, Maybe (Map v Int)) #-}
-doGLP :: (Ord v, Real c) => GLPOpts -> LP v c -> GLPK (ReturnCode, Maybe (Map v Int))
-doGLP SimplexOpts{..} lp = do
-	vars <- writeProblem lp
-	success <- solveSimplex msgLev tmLim presolve
-	bad <- getBadRay
-	maybe (return (success, guard (gaveAnswer success) >> Just vars)) (fail . show) bad
-doGLP MipOpts{..} lp = do
-	vars <- writeProblem lp
-	success <- mipSolve msgLev brTech btTech ppTech fpHeur cuts mipGap tmLim presolve
-	bad <- getBadRay
-	return (success, guard (gaveAnswer success) >> Just vars)
diff --git a/Data/LinearProgram/GLPK/Types.hs b/Data/LinearProgram/GLPK/Types.hs
deleted file mode 100644
--- a/Data/LinearProgram/GLPK/Types.hs
+++ /dev/null
@@ -1,51 +0,0 @@
-{-# LANGUAGE EmptyDataDecls, ForeignFunctionInterface #-}
-
-module Data.LinearProgram.GLPK.Types where
-
-import Control.Monad.Trans (MonadIO (..))
-import Control.Monad (ap)
-
-import Foreign.Ptr
-import Foreign.ForeignPtr
-
-foreign import ccall unsafe "c_glp_create_prob" glpCreateProb :: IO (Ptr GlpProb)
-foreign import ccall unsafe "&c_glp_delete_prob" glpDelProb :: FunPtr (Ptr GlpProb -> IO ())
-
-data GlpProb
-
-data ReturnCode = Success | InvalidBasis | SingularMatrix | IllConditionedMatrix | 
-        InvalidBounds | SolverFailed | ObjLowerLimReached | ObjUpperLimReached | 
-        IterLimReached | TimeLimReached | NoPrimalFeasible | NoDualFeasible | RootLPOptMissing |
-        SearchTerminated | MipGapTolReached | NoPrimDualFeasSolution | NoConvergence |
-        NumericalInstability | InvalidData | ResultOutOfRange deriving (Eq, Show, Enum)
-
-gaveAnswer :: ReturnCode -> Bool
-gaveAnswer = flip elem [Success, IterLimReached, TimeLimReached, SearchTerminated, MipGapTolReached]
-
-newtype GLPK a = GLP {execGLPK :: Ptr GlpProb -> IO a}
-
-runGLPK :: GLPK a -> IO a
-runGLPK m = do  lp <- newForeignPtr glpDelProb =<< glpCreateProb
-                withForeignPtr lp (execGLPK m)
-
-instance Monad GLPK where
-        {-# INLINE return #-}
-        {-# INLINE (>>=) #-}
-        return x = GLP $ \ _ -> return x
-        m >>= k = GLP $ \ lp -> do      x <- execGLPK m lp
-                                        execGLPK (k x) lp
-instance Functor GLPK where
-  fmap f (GLP k) = GLP $ \p -> fmap f (k p)
-
-instance Applicative GLPK where
-  pure = return
-  (<*>) = ap
-
-instance MonadIO GLPK where
-        liftIO m = GLP (const m)
-
-data MsgLev = MsgOff | MsgErr | MsgOn | MsgAll deriving (Eq, Enum, Read, Show)
-data BranchingTechnique = FirstFrac | LastFrac | MostFrac | DrTom | HybridP deriving (Eq, Enum, Read, Show)
-data BacktrackTechnique = DepthFirst | BreadthFirst | LocBound | ProjHeur deriving (Eq, Enum, Read, Show)
-data Preprocessing = NoPre | RootPre | AllPre deriving (Eq, Enum, Read, Show)
-data Cuts = GMI | MIR | Cov | Clq deriving (Eq, Enum, Read, Show)
diff --git a/Data/LinearProgram/LinExpr.hs b/Data/LinearProgram/LinExpr.hs
deleted file mode 100644
--- a/Data/LinearProgram/LinExpr.hs
+++ /dev/null
@@ -1,58 +0,0 @@
-{-# 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
deleted file mode 100644
--- a/Data/LinearProgram/Spec.hs
+++ /dev/null
@@ -1,151 +0,0 @@
-{-# LANGUAGE TupleSections, RecordWildCards, DeriveFunctor #-}
-module Data.LinearProgram.Spec (Constraint(..), VarTypes, ObjectiveFunc, VarBounds, LP(..),
-        mapVars, mapVals, allVars) where
-
-import Control.DeepSeq
-import Control.Monad
-
-import Data.Char (isSpace)
-import Data.Map hiding (map, foldl)
-
-import Text.ParserCombinators.ReadP
-
-import Data.Algebra
-import Data.LinearProgram.Types
-
--- | Representation of a linear constraint on the variables, possibly labeled.
--- The function may be bounded both above and below.
-data Constraint v c = Constr (Maybe String)
-                        (LinFunc v c)
-                        (Bounds c) deriving (Functor)
--- | A mapping from variables to their types.  Variables not mentioned are assumed to be continuous,
-type VarTypes v = Map v VarKind
--- | An objective function for a linear program.
-type ObjectiveFunc = LinFunc
--- | A mapping from variables to their boundaries.  Variables not mentioned are assumed to be free.
-type VarBounds v c = Map v (Bounds c)
-
--- | The specification of a linear programming problem with variables in @v@ and coefficients/constants in @c@.
---   Note: the 'Read' and 'Show' implementations do not correspond to any particular linear program specification format.
-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"
-        Equ x   -> expr ++ " = " ++ show x
-        LBound x -> expr ++ " >= " ++ show x
-        UBound x -> expr ++ " <= " ++ show x
-        Bound l u -> show l ++ " <= " ++ expr ++ " <= " ++ show u
-
-showFunc :: (Show v, Num c, Ord c, Show c) => LinFunc v c -> String
-showFunc func = case assocs func of
-        []      -> "0"
-        ((v,c):vcs) ->
-                show c ++ " " ++ map replaceSpace (show v) ++ 
-                        concatMap showTerm vcs
-        where   showTerm (v, c) = case compare c 0 of
-                        EQ      -> ""
-                        GT      -> " + " ++ show c ++ " " ++ show v
-                        LT      -> " - " ++ show (negate c) ++ " " ++ show v
-                
-replaceSpace :: Char -> Char
-replaceSpace c
-        | isSpace c     = '_'
-        | otherwise     = c
-
-instance (Show v, Num c, Ord c, Show c) => Show (Constraint v c) where
-        show (Constr lab func bds) = maybe "" (++ ": ") lab ++
-                showBds (showFunc func) bds
-
-instance (Read v, Ord v, Read c, Ord c, Num c) => Read (Constraint v c) where
-        readsPrec _= readP_to_S $ liftM toConstr (lab <++ nolab) where
-                toConstr (l, f, bds) = Constr l (fromList f) bds
-                lab = do        skipSpaces
-                                label <- manyTill get (skipSpaces >> char ':')
-                                (_, f, bds) <- nolab
-                                return (Just label, f, bds)
-                nolab = liftM (\ (f, bds) -> (Nothing, f, bds)) $ readBds readConst readFunc
-                readFunc = (do  c <- readCoef readConst
-                                v <- readVar
-                                liftM ((v, c):) readFunc) <++ return []
-                readConst = readS_to_P reads
-                readVar = readS_to_P reads
-
-readCoef :: Num c => ReadP c -> ReadP c
-readCoef readC = between skipSpaces skipSpaces $ 
-        (do     char '+'
-                skipSpaces
-                readC') <++
-        (do     char '-'
-                skipSpaces
-                negate <$> readC') <++ readC'
-        where   readC' = readC <++ return 1
-
-optMaybe :: ReadP a -> ReadP (Maybe a)
-optMaybe p = fmap Just p <++ return Nothing
-
-readBds :: Ord c => ReadP c -> ReadP a -> ReadP (a, Bounds c)
-readBds cst expr = do
-        left <- optMaybe (do    lb <- cst
-                                skipSpaces
-                                rel <- readRelation
-                                return (lb, rel))
-        skipSpaces
-        f <- expr
-        skipSpaces
-        right <- optMaybe (do   rel <- readRelation
-                                skipSpaces
-                                ub <- cst
-                                return (ub, revOrd rel))
-        return (f, getBd left `mappend` getBd right)
-        where   revOrd :: Ordering -> Ordering
-                revOrd GT = LT
-                revOrd LT = GT
-                revOrd EQ = EQ
-                getBd :: Maybe (c, Ordering) -> Bounds c
-                getBd Nothing = Free
-                getBd (Just (x, cmp)) = case cmp of
-                        EQ      -> Equ x
-                        GT      -> LBound x
-                        LT      -> UBound x
-                readRelation = choice [char '<' >> optional (char '=') >> return LT,
-                        char '=' >> return EQ,
-                        char '>' >> optional (char '=') >> return GT]
-
-{-# SPECIALIZE mapVars :: Ord v' => (v -> v') -> LP v Double -> LP v' Double #-}
--- | Applies the specified function to the variables in the linear program.
--- If multiple variables in the original program are mapped to the same variable in the new program,
--- in general, we set those variables to all be equal, as follows.
--- 
--- * In linear functions, including the objective function and the constraints,
---      coefficients will be added together.  For instance, if @v1,v2@ are mapped to the same
---      variable @v'@, then a linear function of the form @c1 *& v1 ^+^ c2 *& v2@ will be mapped to
---      @(c1 ^+^ c2) *& v'@.
---
--- * In variable bounds, bounds will be combined.  An error will be thrown if the bounds
---      are mutually contradictory.
--- 
--- * In variable kinds, the most restrictive kind will be retained.
-mapVars :: (Ord v', Ord c, Group c) => (v -> v') -> LP v c -> LP v' c
-mapVars f LP{..} =  
-        LP{objective = mapKeysWith (^+^) f objective, 
-                constraints = [Constr lab (mapKeysWith (^+^) f func) bd | Constr lab func bd <- constraints],
-                varBounds = mapKeysWith mappend f varBounds,
-                varTypes = mapKeysWith mappend f varTypes, ..}
-
--- | Applies the specified function to the constants in the linear program.  This is only safe
--- for a monotonic function.
-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 (LP v c) where
-        rnf LP{..} = direction `deepseq` objective `deepseq` constraints `deepseq`
-                varBounds `deepseq` rnf varTypes
diff --git a/Data/LinearProgram/Types.hs b/Data/LinearProgram/Types.hs
deleted file mode 100644
--- a/Data/LinearProgram/Types.hs
+++ /dev/null
@@ -1,72 +0,0 @@
-{-# LANGUAGE DeriveFunctor, DeriveGeneric #-}
-module Data.LinearProgram.Types (VarKind(..), Direction(..), Bounds(..)) where
-
-import Control.DeepSeq
-import Data.Monoid
-import GHC.Generics
-
-data VarKind = ContVar | IntVar | BinVar deriving (Eq, Ord, Enum, Show, Read, Generic)
-
--- instance NFData VarKind
-
-instance Monoid VarKind where
-        mempty = ContVar
-        mappend = max
-
-data Direction = Min | Max deriving (Eq, Ord, Enum, Show, Read, Generic)
-
--- instance NFData Direction
-
-data Bounds a =
-        Free | LBound !a | UBound !a | Equ !a | Bound !a !a deriving (Eq, Show, Read, Functor)
-
-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 (Bounds a)
-
--- Bounds form a monoid under intersection.
-instance Ord a => Monoid (Bounds a) where
-        mempty = Free
-        Free `mappend` bd = bd
-        bd `mappend` Free = bd
-        Equ a `mappend` Equ b
-                | a == b        = Equ a
-        Equ a `mappend` UBound b
-                | a <= b        = Equ a
-        Equ a `mappend` LBound b
-                | a >= b        = Equ a
-        Equ a `mappend` Bound l u
-                | a >= l && a <= u
-                                = Equ a
-        Equ _ `mappend` _ = infeasible
-        UBound b `mappend` Equ a
-                | a <= b        = Equ a
-        LBound b `mappend` Equ a
-                | a >= b        = Equ a
-        Bound l u `mappend` Equ a
-                | a >= l && a <= u
-                                = Equ a
-        _ `mappend` Equ _ = infeasible
-        LBound a `mappend` LBound b = LBound (max a b)
-        LBound l `mappend` UBound u = bound l u
-        UBound u `mappend` LBound l = bound l u
-        LBound a `mappend` Bound l u = bound (max a l) u
-        Bound l u `mappend` LBound a = bound (max a l) u
-        UBound a `mappend` UBound b = UBound (min a b)
-        UBound a `mappend` Bound l u = bound l (min a u)
-        Bound l u `mappend` UBound a = bound l (min a u)
-        Bound l u `mappend` Bound l' u' = bound (max l l') (min u u')
-
-infeasible :: Bounds a
-infeasible = error "Mutually contradictory constraints found."
-
-bound :: Ord a => a -> a -> Bounds a
-bound l u       | l <= u        = Bound l u
-                | otherwise     = infeasible
diff --git a/Setup.hs b/Setup.hs
deleted file mode 100644
--- a/Setup.hs
+++ /dev/null
@@ -1,2 +0,0 @@
-import Distribution.Simple
-main = defaultMain
diff --git a/examples/example1.hs b/examples/example1.hs
--- a/examples/example1.hs
+++ b/examples/example1.hs
@@ -1,21 +1,29 @@
+import Prelude hiding (Num(..))
 
-import Data.LinearProgram.LPMonad
+import Algebra.Classes
+import Control.Monad.LPMonad
+import Data.LinearProgram.Common
 import Data.LinearProgram
 import Data.LinearProgram.GLPK
+import qualified Data.Map as M
+import Data.LinearProgram.LinExpr
 
 objFun :: LinFunc String Int
 objFun = linCombination [(10, "x1"), (6, "x2"), (4, "x3")]
 
+n *& v = linCombination [(n,v)]
+
 lp :: LP String Int
-lp = execLPM $ do	setDirection Max
-			setObjective objFun
-			leqTo (varSum ["x1", "x2", "x3"]) 100
-			leqTo (10 *^ var "x1" ^+^ 4 *& "x2" ^+^ 5 *^ var "x3") 600
-			leqTo (linCombination [(2, "x1"), (2, "x2"), (6, "x3")]) 300
-			varGeq "x1" 0
-			varBds "x2" 0 50
-			varGeq "x3" 0
-			setVarKind "x1" IntVar
-			setVarKind "x2" ContVar
+lp = execLPM $ do
+  setDirection Max
+  setObjective objFun
+  leqTo (Data.LinearProgram.sum $ map (1 *&) ["x1", "x2", "x3"]) 100
+  leqTo (10 *& "x1" + 4 *& "x2" + 5 *& "x3") 600
+  leqTo (linCombination [(2, "x1"), (2, "x2"), (6, "x3")]) 300
+  varGeq "x1" 0
+  varBds "x2" 0 50
+  varGeq "x3" 0
+  setVarKind "x1" IntVar
+  setVarKind "x2" ContVar
 
 main = print =<< glpSolveVars mipDefaults lp
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.3.5
+Version:        0.8
 Author:         Louis Wasserman
 License:        BSD3
 License-file:   LICENSE
@@ -17,24 +17,28 @@
     of options available.
 
 Category:      Math
-cabal-version: >= 1.6
-build-type:     Simple
-
-extra-source-files: examples/example1.hs
+cabal-version: 1.12
+build-type:    Simple
 
 source-repository head
   type: git
   location: https://github.com/jyp/glpk-hs
 
+executable glpk-hs-example
+  main-is:          examples/example1.hs
+  build-depends:    base >= 4 && < 5, array, containers, mtl, deepseq, gasp, glpk-hs
+  ghc-options:      -O2 -Wall
+  default-language: Haskell2010
+
 library
-  Build-Depends:    base >= 4 && < 5, array, containers, mtl, deepseq
+  default-language: Haskell2010
+  Build-Depends:    base >= 4 && < 5, array, containers, mtl, deepseq, gasp >= 1.2
   Exposed-modules:  Data.LinearProgram,
                     Data.LinearProgram.Common,
                     Data.LinearProgram.LinExpr,
                     Data.LinearProgram.GLPK,
                     Data.LinearProgram.GLPK.Solver,
                     Data.LinearProgram.GLPK.IO,
-                    Data.Algebra,
                     Control.Monad.LPMonad,
                     Control.Monad.LPMonad.Supply,
                     Control.Monad.LPMonad.Supply.Class
@@ -44,10 +48,15 @@
                     Data.LinearProgram.GLPK.IO.Internal,
                     Control.Monad.LPMonad.Internal,
                     Data.LinearProgram.Spec,
-                    Data.LinearProgram.Types,
-                    Data.Algebra.Group,
-                    Data.Algebra.Ring,
-                    Data.Algebra.Module,
-                    Data.Algebra.Field
+                    Data.LinearProgram.Types
+  hs-source-dirs:   src
   c-sources:        glpk/glpk.c
   extra-libraries:  glpk
+  if os(OSX)
+      extra-lib-dirs: /usr/lib
+      extra-lib-dirs: /opt/local/lib/
+      include-dirs: /opt/local/include/
+      extra-lib-dirs: /usr/local/lib/
+      include-dirs: /usr/local/include/
+      if arch(i386)
+          cc-options: -arch i386
diff --git a/src/Control/Monad/LPMonad.hs b/src/Control/Monad/LPMonad.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/LPMonad.hs
@@ -0,0 +1,98 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+-- | A collection of operations that can be used to specify linear programming in a
+-- simple, monadic way.  It is not too difficult to construct 'LP' values explicitly,
+-- but this module may help simplify and modularize the construction of the linear program,
+-- for example separating different families of constraints in the problem specification.
+-- 
+-- Many of these functions should be executed in either the @'LPM' v c@ or the @'LPT' v c 'IO'@ monad.
+-- If you wish to generate new variables on an ad-hoc basis, rather than supplying your own variable type, use the
+-- 'VSupply' or 'VSupplyT' monads in your transformer stack, as in @'LPT' 'Var' c 'VSupply'@ or
+-- @'LPT' 'Var' c ('VSupplyT' 'IO')@.  To generate new variables, use 'supplyNew' or 'supplyN'.
+module Control.Monad.LPMonad (
+	module Control.Monad.LPMonad.Internal,
+	-- * Generation of new variables
+	module Control.Monad.LPMonad.Supply,
+	-- * Solvers
+	quickSolveMIP,
+	quickSolveLP,
+	glpSolve,
+	quickSolveMIP',
+	quickSolveLP',
+	glpSolve',
+	-- * File I/O
+	writeLPToFile,
+	readLPFromFile,
+	readLPFromFile') where
+
+import Control.Monad ((<=<))
+import Control.Monad.State.Class (MonadState(..))
+import Control.Monad.Trans (MonadIO (..))
+
+import Data.Map (Map)
+
+import Data.LinearProgram.Common
+import Control.Monad.LPMonad.Internal
+import Control.Monad.LPMonad.Supply
+
+import Data.LinearProgram.GLPK.Solver
+import Data.LinearProgram.GLPK.IO
+
+{-# SPECIALIZE quickSolveLP :: (Ord v, Real c) => 
+	LPT v c IO (ReturnCode, Maybe (Double, Map v Double)) #-}
+{-# SPECIALIZE quickSolveMIP :: (Ord v, Real c) => 
+	LPT v c IO (ReturnCode, Maybe (Double, Map v Double)) #-}
+-- | Solves the linear program with the default settings in GLPK.  Returns the return code,
+-- and if the solver was successful, the objective function value and the settings of each variable.
+quickSolveLP, quickSolveMIP :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => 
+	m (ReturnCode, Maybe (Double, Map v Double))
+quickSolveLP = glpSolve simplexDefaults
+quickSolveMIP = glpSolve mipDefaults
+
+{-# SPECIALIZE glpSolve :: (Ord v, Real c) => GLPOpts -> LPT v c IO (ReturnCode, Maybe (Double, Map v Double)) #-}
+-- | Solves the linear program with the specified options in GLPK.  Returns the return code,
+-- and if the solver was successful, the objective function value and the settings of each variable.
+glpSolve :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => GLPOpts -> m (ReturnCode, Maybe (Double, Map v Double))
+glpSolve opts = get >>= liftIO . glpSolveVars opts
+
+{-# SPECIALIZE quickSolveLP' :: (Ord v, Real c) => LPT v c IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c])) #-}
+{-# SPECIALIZE quickSolveMIP' :: (Ord v, Real c) => LPT v c IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c])) #-}
+-- | Solves the linear program with the default settings in GLPK.  Returns the return code,
+-- and if the solver was successful, the objective function value, the settings of each variable, and the
+-- value of each constraint/row.
+quickSolveLP', quickSolveMIP' :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => 
+	m (ReturnCode, Maybe (Double, Map v Double, [RowValue v c]))
+quickSolveLP' = glpSolve' simplexDefaults
+quickSolveMIP' = glpSolve' mipDefaults
+
+{-# SPECIALIZE glpSolve' :: (Ord v, Real c) => GLPOpts -> LPT v c IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c])) #-}
+-- | Solves the linear program with the specified options in GLPK.  Returns the return code,
+-- and if the solver was successful, the objective function value, the settings of each variable, and
+-- the value of each constraint/row.
+glpSolve' :: (Ord v, Real c, MonadState (LP v c) m, MonadIO m) => 
+	GLPOpts -> m (ReturnCode, Maybe (Double, Map v Double, [RowValue v c]))
+glpSolve' opts = get >>= liftIO . glpSolveAll opts
+
+{-# SPECIALIZE writeLPToFile :: (Ord v, Show v, Real c) => FilePath -> LPT v c IO () #-}
+-- | Writes the current linear program to the specified file in CPLEX LP format. 
+-- (This is a binding to GLPK, not a Haskell implementation of CPLEX.)
+writeLPToFile :: (Ord v, Show v, Real c, MonadState (LP v c) m, MonadIO m) =>
+	FilePath -> m ()
+writeLPToFile file = get >>= liftIO . writeLP file 
+
+{-# SPECIALIZE readLPFromFile :: (Ord v, Read v, Fractional c) => FilePath -> LPT v c IO () #-}
+-- | Reads a linear program from the specified file in CPLEX LP format, overwriting
+-- the current linear program.  Uses 'read' and 'realToFrac' to translate to the specified type.
+-- Warning: this may not work on all files written using 'writeLPToFile', since variable names
+-- may be changed.
+-- (This is a binding to GLPK, not a Haskell implementation of CPLEX.)
+readLPFromFile :: (Ord v, Read v, Fractional c, MonadState (LP v c) m, MonadIO m) =>
+	FilePath -> m ()
+readLPFromFile = put <=< liftIO . readLP
+
+{-# SPECIALIZE readLPFromFile :: FilePath -> LPT String Double IO () #-}
+-- | Reads a linear program from the specified file in CPLEX LP format, overwriting
+-- the current linear program.  (This is a binding to GLPK, not a Haskell implementation of CPLEX.)
+readLPFromFile' :: (MonadState (LP String Double) m, MonadIO m) =>
+	FilePath -> m ()
+readLPFromFile' = put <=< liftIO . readLP'
diff --git a/src/Control/Monad/LPMonad/Internal.hs b/src/Control/Monad/LPMonad/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/LPMonad/Internal.hs
@@ -0,0 +1,248 @@
+{-# LANGUAGE BangPatterns, FlexibleContexts, RecordWildCards #-}
+
+module Control.Monad.LPMonad.Internal (
+--         module Data.LinearProgram.Common,
+        -- * Monad definitions
+        LPM,
+        LPT,
+        runLPM,
+        runLPT,
+        execLPM,
+        execLPT,
+        evalLPM,
+        evalLPT,
+        -- * Constructing the LP
+        -- ** Objective configuration
+        setDirection,
+        setObjective,
+        addObjective,
+        addWeightedObjective,
+        -- ** Two-function constraints
+        leq,
+        equal,
+        geq,
+        leq',
+        equal',
+        geq',
+        -- ** One-function constraints
+        leqTo,
+        equalTo,
+        geqTo,
+        constrain,
+        leqTo',
+        equalTo',
+        geqTo',
+        constrain',
+        -- ** Variable constraints
+        varLeq,
+        varEq,
+        varGeq,
+        varBds,
+        setVarBounds,
+        setVarKind,
+--         newVariables,
+--         newVariables'
+        ) where
+
+import Prelude hiding ((-),(+))
+import Control.Monad.State.Strict
+import Control.Monad.Identity
+
+import Data.Map
+
+import Data.LinearProgram.Common
+
+-- | A simple monad for constructing linear programs.  This library is intended to be able to link to
+-- a variety of different linear programming implementations.
+type LPM v c = LPT v c Identity
+
+-- | A simple monad transformer for constructing linear programs in an arbitrary monad.
+type LPT v c = StateT (LP v c)
+
+runLPM :: (Ord v, Group c) => LPM v c a -> (a, LP v c)
+runLPM = runIdentity . runLPT
+
+runLPT :: (Ord v, Group c) => LPT v c m a -> m (a, LP v c)
+runLPT m = runStateT m (LP Max zero [] mempty mempty)
+
+-- | Constructs a linear programming problem.
+execLPM :: (Ord v, Group c) => LPM v c a -> LP v c
+execLPM = runIdentity . execLPT
+
+-- | Constructs a linear programming problem in the specified monad.
+execLPT :: (Ord v, Group c, Monad m) => LPT v c m a -> m (LP v c)
+execLPT = liftM snd . runLPT
+
+-- | Runs the specified operation in the linear programming monad.
+evalLPM :: (Ord v, Group c) => LPM v c a -> a
+evalLPM = runIdentity . evalLPT
+
+-- | Runs the specified operation in the linear programming monad transformer.
+evalLPT :: (Ord v, Group c, Monad m) => LPT v c m a -> m a
+evalLPT = liftM fst . runLPT
+
+-- | Sets the optimization direction of the linear program: maximization or minimization.
+{-# SPECIALIZE setDirection :: Direction -> LPM v c (), Monad m => Direction -> LPT v c m () #-}
+setDirection :: (MonadState (LP v c) m) => Direction -> m ()
+setDirection dir = modify (\ lp -> lp{direction = dir})
+
+{-# SPECIALIZE equal :: (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 () #-}
+{-# SPECIALIZE leq :: (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 () #-}
+{-# 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.  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
+geq = flip leq
+
+{-# SPECIALIZE equal' :: (Ord v, Group c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
+        (Ord v, Group c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+{-# SPECIALIZE geq' :: (Ord v, Group c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
+        (Ord v, Group c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+{-# SPECIALIZE leq' :: (Ord v, Group c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
+        (Ord v, Group c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+-- | Specifies the relationship between two functions in the variables, with a label on the constraint.
+equal', leq', geq' :: (Ord v, Group c, MonadState (LP v c) m) => String -> LinFunc v c -> LinFunc v c -> m ()
+equal' lab f g = equalTo' lab (f - g) zero
+leq' lab f g = leqTo' lab (f - g) zero
+geq' = flip . leq'
+
+{-# SPECIALIZE equalTo :: LinFunc v c -> c -> LPM v c (), Monad m => LinFunc v c -> c -> LPT v c m () #-}
+{-# SPECIALIZE geqTo :: LinFunc v c -> c -> LPM v c (), Monad m => LinFunc v c -> c -> LPT v c m () #-}
+{-# SPECIALIZE leqTo :: LinFunc v c -> c -> LPM v c (), Monad m => LinFunc v c -> c -> LPT v c m () #-}
+-- | Sets a constraint on a linear function in the variables.
+equalTo, leqTo, geqTo :: MonadState (LP v c) m => LinFunc v c -> c -> m ()
+equalTo f v = constrain f (Equ v)
+leqTo f v = constrain f (UBound v)
+geqTo f v = constrain f (LBound v)
+
+{-# SPECIALIZE equalTo' :: String -> LinFunc v c -> c -> LPM v c (),
+        Monad m => String -> LinFunc v c -> c -> LPT v c m () #-}
+{-# SPECIALIZE geqTo' :: String -> LinFunc v c -> c -> LPM v c (),
+        Monad m => String -> LinFunc v c -> c -> LPT v c m () #-}
+{-# SPECIALIZE leqTo' :: String -> LinFunc v c -> c -> LPM v c (),
+        Monad m => String -> LinFunc v c -> c -> LPT v c m () #-}
+-- | Sets a labeled constraint on a linear function in the variables.
+equalTo', leqTo', geqTo' :: MonadState (LP v c) m => String -> LinFunc v c -> c -> m ()
+equalTo' lab f v = constrain' lab f (Equ v)
+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 = () <$ 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 () #-}
+{-# SPECIALIZE varLeq :: (Ord v, Ord c) => v -> c -> LPM v c (),
+        (Ord v, Ord c, Monad m) => v -> c -> LPT v c m () #-}
+{-# SPECIALIZE varGeq :: (Ord v, Ord c) => v -> c -> LPM v c (),
+        (Ord v, Ord c, Monad m) => v -> c -> LPT v c m () #-}
+-- | Sets a constraint on the value of a variable.  If you constrain a variable more than once,
+-- the constraints will be combined.  If the constraints are mutually contradictory,
+-- an error will be generated.  This is more efficient than adding an equivalent function constraint.
+varEq, varLeq, varGeq :: (Ord v, Ord c, MonadState (LP v c) m) => v -> c -> m ()
+varEq v c = setVarBounds v (Equ c)
+varLeq v c = setVarBounds v (UBound c)
+varGeq v c = setVarBounds v (LBound c)
+
+{-# SPECIALIZE varBds :: (Ord v, Ord c) => v -> c -> c -> LPM v c (),
+        (Ord v, Ord c, Monad m) => v -> c -> c -> LPT v c m () #-}
+-- | Bounds the value of a variable on both sides.  If you constrain a variable more than once,
+-- the constraints will be combined.  If the constraints are mutually contradictory,
+-- an error will be generated.  This is more efficient than adding an equivalent function constraint.
+varBds :: (Ord v, Ord c, MonadState (LP v c) m) => v -> c -> c -> m ()
+varBds v l u = setVarBounds v (Bound l u)
+
+{-# SPECIALIZE constrain :: LinFunc v c -> Bounds c -> LPM v c (),
+        Monad m => LinFunc v c -> Bounds c -> LPT v c m () #-}
+-- | The most general form of an unlabeled constraint.
+constrain :: MonadState (LP v c) m => LinFunc v c -> Bounds c -> m ()
+constrain f bds = modify addConstr where
+        addConstr lp@LP{..}
+                = lp{constraints = Constr Nothing f bds:constraints}
+
+{-# SPECIALIZE constrain' :: String -> LinFunc v c -> Bounds c -> LPM v c (),
+        Monad m => String -> LinFunc v c -> Bounds c -> LPT v c m () #-}
+-- | The most general form of a labeled constraint.
+constrain' :: MonadState (LP v c) m => String -> LinFunc v c -> Bounds c -> m ()
+constrain' lab f bds = modify addConstr where
+        addConstr lp@LP{..}
+                = lp{constraints = Constr (Just lab) f bds:constraints}
+
+{-# SPECIALIZE setObjective :: LinFunc v c -> LPM v c (),
+        Monad m => LinFunc v c -> LPT v c m () #-}
+-- | Sets the objective function, overwriting the previous objective function.
+setObjective :: MonadState (LP v c) m => LinFunc v c -> m ()
+setObjective obj = modify setObj where
+        setObj lp = lp{objective = obj}
+
+{-# SPECIALIZE addObjective :: (Ord v, Group c) => LinFunc v c -> LPM v c (),
+        (Ord v, Group c, Monad m) => LinFunc v c -> LPT v c m () #-}
+-- | Adds this function to the objective function.
+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, Ring c) => c -> LinFunc v c -> LPM v c (),
+        (Ord v, Ring c, Monad m) => c -> LinFunc v c -> LPT v c m () #-}
+-- | Adds this function to the objective function, with the specified weight.  Equivalent to
+-- @'addObjective' (wt '*^' obj)@.
+addWeightedObjective :: (Ord v, Ring c, MonadState (LP v c) m) =>
+                        c -> LinFunc v c -> m ()
+addWeightedObjective wt obj = addObjective (wt *^ obj)
+
+{-# SPECIALIZE setVarBounds :: (Ord v, Ord c) => v -> Bounds c -> LPM v c (),
+        (Ord v, Ord c, Monad m) => v -> Bounds c -> LPT v c m () #-}
+-- | The most general way to set constraints on a variable.
+-- If you constrain a variable more than once, the constraints will be combined.
+-- If you combine mutually contradictory constraints, an error will be generated.
+-- This is more efficient than creating an equivalent function constraint.
+setVarBounds :: (Ord v, Ord c, MonadState (LP v c) m) => v -> Bounds c -> m ()
+setVarBounds var bds = modify addBds where
+        addBds lp@LP{..} = lp{varBounds = insertWith mappend var bds varBounds}
+
+{-# SPECIALIZE setVarKind :: Ord v => v -> VarKind -> LPM v c (),
+        (Ord v, Monad m) => v -> VarKind -> LPT v c m () #-}
+-- | Sets the kind ('type') of a variable.  See 'VarKind'.
+setVarKind :: (Ord v, MonadState (LP v c) m) => v -> VarKind -> m ()
+setVarKind v k = modify setK where
+        setK lp@LP{..} = lp{varTypes = insertWith mappend v k varTypes}
diff --git a/src/Control/Monad/LPMonad/Supply.hs b/src/Control/Monad/LPMonad/Supply.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/LPMonad/Supply.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}
+
+module Control.Monad.LPMonad.Supply (module Control.Monad.LPMonad.Supply.Class, Var(..), VSupply, VSupplyT, runVSupply, runVSupplyT) where
+
+import Control.Monad.Identity
+import Control.Monad.Trans
+import Control.Monad.State.Strict
+import Control.Monad.RWS.Class
+import Control.Monad.Cont.Class
+import Control.Monad.Error.Class
+import Control.Applicative
+import Control.Monad.LPMonad.Supply.Class
+
+-- | A type suitable for use as a linear program variable.
+newtype Var = Var {varId :: Int} deriving (Eq, Ord, Enum)
+
+-- | A monad capable of supplying unique variables.
+type VSupply = VSupplyT Identity
+
+runVSupply :: VSupply a -> a
+runVSupply = runIdentity . runVSupplyT
+
+-- | A monad transformer capable of supplying unique variables.
+newtype VSupplyT m a = VSupplyT (StateT Var m a) deriving (Functor, Applicative, Monad, Alternative, MonadPlus, MonadTrans, MonadReader r, MonadWriter w, MonadCont,
+        MonadIO, MonadFix, MonadError e)
+
+runVSupplyT :: Monad m => VSupplyT m a -> m a
+runVSupplyT (VSupplyT m) = evalStateT m (Var 0)
+
+instance Show Var where
+        show (Var x) = "x_" ++ show x
+
+instance Read Var where
+        readsPrec _ ('x':'_':xs) = [(Var x, s') | (x, s') <- reads xs]
+        readsPrec _ _ = []
+
+instance MonadState s m => MonadState s (VSupplyT m) where
+        get = lift get
+        put = lift . put
+
+instance Monad m => MonadSupply Var (VSupplyT m) where
+        {-# SPECIALIZE instance MonadSupply Var VSupply #-}
+        supplyNew = VSupplyT $ StateT $ \ v -> return (v, succ v)
+        supplyN n = VSupplyT $ StateT $ \ (Var x) -> return (map Var [x..x+n-1], Var (x + n))
diff --git a/src/Control/Monad/LPMonad/Supply/Class.hs b/src/Control/Monad/LPMonad/Supply/Class.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/LPMonad/Supply/Class.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE UndecidableInstances, FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies #-}
+module Control.Monad.LPMonad.Supply.Class where
+
+import Control.Monad
+
+import Control.Monad.State.Strict
+import Control.Monad.Reader
+import Control.Monad.Error
+import qualified Control.Monad.Writer.Lazy as WL
+import qualified Control.Monad.Writer.Strict as WS
+import qualified Control.Monad.State.Lazy as SL
+import Control.Monad.Cont
+
+import Data.Monoid
+
+-- | A class implemented by monads that can supply values of type @s@.  Minimal implementation: 'supplyNew' or 'supplyN'.
+class Monad m => MonadSupply s m | m -> s where
+	-- | Supply a new value of type @s@.
+	supplyNew :: m s
+	-- | Supply @n@ values of type @s@.
+	supplyN :: Int -> m [s]
+	
+	supplyNew = liftM head (supplyN 1)
+	supplyN n = replicateM n supplyNew
+
+instance MonadSupply x m => MonadSupply x (StateT s m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
+
+instance MonadSupply x m => MonadSupply x (ReaderT r m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
+
+instance (Error e, MonadSupply x m) => MonadSupply x (ErrorT e m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
+
+instance (MonadSupply x m, Monoid w) => MonadSupply x (WL.WriterT w m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
+
+instance (MonadSupply x m, Monoid w) => MonadSupply x (WS.WriterT w m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
+
+instance MonadSupply x m => MonadSupply x (ContT r m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
+
+instance MonadSupply x m => MonadSupply x (SL.StateT s m) where
+	supplyNew = lift supplyNew
+	supplyN = lift . supplyN
diff --git a/src/Data/LinearProgram.hs b/src/Data/LinearProgram.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram.hs
@@ -0,0 +1,8 @@
+module Data.LinearProgram (
+	module Data.LinearProgram.Common,
+	module Data.LinearProgram.GLPK,
+	module Control.Monad.LPMonad) where
+
+import Data.LinearProgram.GLPK
+import Data.LinearProgram.Common
+import Control.Monad.LPMonad
diff --git a/src/Data/LinearProgram/Common.hs b/src/Data/LinearProgram/Common.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/Common.hs
@@ -0,0 +1,19 @@
+-- | Contains sufficient tools to represent linear programming problems in Haskell.  In the future, if linkings to other
+-- linear programming libraries are made, this will be common to them all.
+module Data.LinearProgram.Common (
+	module Data.LinearProgram.Spec,
+	module Algebra.Classes,
+	module Data.LinearProgram.Types) where
+
+import Data.LinearProgram.Spec
+import Algebra.Classes
+import Data.LinearProgram.Types
+
+import Data.Map
+import GHC.Exts (build)
+
+{-# RULES
+	"assocs" assocs = \ m -> build (\ c n -> foldrWithKey (curry c) n m);
+	"elems" elems = \ m -> build (\ c n -> foldrWithKey (const c) n m);
+	"keys" keys = \ m -> build (\ c n -> foldrWithKey (\ k _ -> c k) n m);
+	#-}
diff --git a/src/Data/LinearProgram/GLPK.hs b/src/Data/LinearProgram/GLPK.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK.hs
@@ -0,0 +1,8 @@
+module Data.LinearProgram.GLPK (
+-- 	module Data.LinearProgram.LPMonad,
+	module Data.LinearProgram.GLPK.Solver,
+	module Data.LinearProgram.GLPK.IO) where
+	
+-- import Data.LinearProgram.LPMonad
+import Data.LinearProgram.GLPK.Solver
+import Data.LinearProgram.GLPK.IO
diff --git a/src/Data/LinearProgram/GLPK/Common.hs b/src/Data/LinearProgram/GLPK/Common.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK/Common.hs
@@ -0,0 +1,13 @@
+module Data.LinearProgram.GLPK.Common (
+	module Data.LinearProgram.GLPK.Internal,
+	module Data.LinearProgram.GLPK.Types,
+	module Foreign.Ptr,
+	module Foreign.C,
+	module Foreign.Marshal.Array) where
+
+import Data.LinearProgram.GLPK.Internal
+import Data.LinearProgram.GLPK.Types
+
+import Foreign.Ptr
+import Foreign.C
+import Foreign.Marshal.Array
diff --git a/src/Data/LinearProgram/GLPK/IO.hs b/src/Data/LinearProgram/GLPK/IO.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK/IO.hs
@@ -0,0 +1,21 @@
+-- | Bindings to the file I/O functions from GLPK, on the CPLEX LP file format.
+module Data.LinearProgram.GLPK.IO where
+
+import Data.LinearProgram.Common
+
+import Data.LinearProgram.GLPK.Common
+import Data.LinearProgram.GLPK.IO.Internal
+
+{-# SPECIALIZE readLP :: (Ord v, Read v) => FilePath -> IO (LP v Double) #-}
+-- | Read a linear program from a file in CPLEX LP format.  Warning: this will not necessarily succeed
+-- on all files generated by 'writeLP', as variable names may be changed.
+readLP :: (Ord v, Read v, Fractional c) => FilePath -> IO (LP v c)
+readLP = fmap (mapVals realToFrac . mapVars read) . readLP'
+
+-- | Read a linear program from a file in CPLEX LP format.
+readLP' :: FilePath -> IO (LP String Double)
+readLP' = runGLPK . readGLPLP
+
+-- | Write a linear program to a file in CPLEX LP format.
+writeLP :: (Ord v, Show v, Real c) => FilePath -> LP v c -> IO ()
+writeLP file = runGLPK . writeGLPLP file
diff --git a/src/Data/LinearProgram/GLPK/IO/Internal.hs b/src/Data/LinearProgram/GLPK/IO/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK/IO/Internal.hs
@@ -0,0 +1,134 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+module Data.LinearProgram.GLPK.IO.Internal (readGLPLP, writeGLPLP) where
+import Prelude hiding ((+))
+import Control.Monad
+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
+import Control.Monad.LPMonad.Internal
+
+foreign import ccall unsafe "c_glp_write_lp" glpWriteLP :: Ptr GlpProb -> CString -> IO ()
+foreign import ccall unsafe "c_glp_read_lp" glpReadLP :: Ptr GlpProb -> CString -> IO ()
+foreign import ccall unsafe "c_glp_set_col_name" glpSetColName :: Ptr GlpProb -> CInt -> CString -> IO ()
+foreign import ccall unsafe "c_glp_set_row_name" glpSetRowName :: Ptr GlpProb -> CInt -> CString -> IO ()
+foreign import ccall unsafe "c_glp_get_obj_dir" glpGetObjDir :: Ptr GlpProb -> IO CInt
+foreign import ccall unsafe "c_glp_get_num_rows" glpGetNumRows :: Ptr GlpProb -> IO CInt
+foreign import ccall unsafe "c_glp_get_num_cols" glpGetNumCols :: Ptr GlpProb -> IO CInt
+foreign import ccall unsafe "c_glp_get_row_name" glpGetRowName :: Ptr GlpProb -> CInt -> IO CString
+foreign import ccall unsafe "c_glp_get_col_name" glpGetColName :: Ptr GlpProb -> CInt -> IO CString
+foreign import ccall unsafe "c_glp_get_col_kind" glpGetColKind :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_get_row_type" glpGetRowType :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_get_col_type" glpGetColType :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_get_row_lb" glpGetRowLb :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_get_col_lb" glpGetColLb :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_get_row_ub" glpGetRowUb :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_get_col_ub" glpGetColUb :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_get_obj_coef" glpGetObjCoef :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_get_mat_row" glpGetMatRow :: Ptr GlpProb -> CInt -> Ptr CInt -> Ptr CDouble -> IO CInt
+
+writeLP :: FilePath -> GLPK ()
+writeLP file = GLP $ withCString file . glpWriteLP
+
+readLP :: FilePath -> GLPK ()
+readLP file = GLP $ withCString file . glpReadLP
+
+getDir :: GLPK Direction
+getDir = liftM (toEnum . subtract 1 . fromIntegral) $ GLP glpGetObjDir
+
+getRowName, getColName :: Int -> GLPK (Maybe String)
+getRowName i = GLP $ peekCAString' <=< flip glpGetRowName (fromIntegral i)
+getColName i = GLP $ peekCAString' <=< flip glpGetColName (fromIntegral i)
+
+peekCAString' :: CString -> IO (Maybe String)
+peekCAString' str
+	| str == nullPtr	= return Nothing
+	| otherwise		= liftM Just $ peekCAString str
+
+getNumRows, getNumCols :: GLPK Int
+getNumRows = liftM fromIntegral $ GLP glpGetNumRows
+getNumCols = liftM fromIntegral $ GLP glpGetNumCols
+
+rowBounds, colBounds :: Int -> GLPK (Bounds Double)
+rowBounds = loadBounds (getCDouble glpGetRowLb) (getCDouble glpGetRowUb) (getCInt glpGetRowType)
+colBounds = loadBounds (getCDouble glpGetColLb) (getCDouble glpGetColUb) (getCInt glpGetColType)
+
+colKind :: Int -> GLPK VarKind
+colKind = liftM (toEnum . subtract 1) . getCInt glpGetColKind
+
+getCInt :: (Ptr GlpProb -> CInt -> IO CInt) -> Int -> GLPK Int
+getCInt f i = GLP $ \ lp -> liftM fromIntegral $ f lp (fromIntegral i)
+
+getCDouble :: (Ptr GlpProb -> CInt -> IO CDouble) -> Int -> GLPK Double
+getCDouble f i = GLP $ \ lp -> liftM realToFrac $ f lp (fromIntegral i)
+
+setRowName :: Int -> String -> GLPK ()
+setRowName i nam = GLP $ withCString nam . flip glpSetRowName (fromIntegral i)
+
+setColName :: Int -> String -> GLPK ()
+setColName i nam = GLP $ withCString nam . flip glpSetColName (fromIntegral i)
+
+loadBounds :: (Int -> GLPK Double) -> (Int -> GLPK Double) ->
+	(Int -> GLPK Int) -> Int -> GLPK (Bounds Double)
+loadBounds lb ub tp i = do
+	typ <- tp i
+	case typ of
+		1	-> return Free
+		2	-> liftM LBound (lb i)
+		3	-> liftM UBound (ub i)
+		4	-> liftM2 Bound (lb i) (ub i)
+		_	-> liftM Equ (lb i)
+
+getObjCoef :: Int -> GLPK Double
+getObjCoef = getCDouble glpGetObjCoef
+
+getRows :: GLPK [(Int, [(Int, Double)])]
+getRows = do	n <- getNumRows
+		m <- getNumCols
+		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 $ mapM (peekElemOff ixs) [1..k]
+			coefsL <- liftIO $ mapM (peekElemOff ixs) [1..k]
+			return (i, zip (map fromIntegral ixsL) (map realToFrac coefsL))
+			| i <- [1..n]]
+
+readGLPLP :: FilePath -> GLPK (LP String Double)
+readGLPLP file = execLPT $ do
+	lift $ readLP file
+	setDirection =<< lift getDir
+	nCols <- lift getNumCols
+	names <- lift $ liftM fromList $ mapM (\ i -> do
+		Just name <- getColName i
+		return (i, name)) [1..nCols]
+	sequence_ [do
+		bds <- lift $ colBounds i
+		kind <- lift $ colKind i
+		setVarBounds name bds
+		setVarKind name kind
+		return (i, name)
+			| (i, name) <- assocs names]
+	rowContents <- lift getRows
+	sequence_ [do
+		bds <- lift $ rowBounds i
+		name <- lift $ getRowName i
+		maybe constrain constrain' name
+			(linCombination [(v, names ! j) | (j, v) <- row]) bds
+			| (i, row) <- rowContents]
+	obj <- lift $ sequence [do
+		c <- getObjCoef i
+		return (name, c) | (i, name) <- assocs names]
+	setObjective (fromList (filter ((/= 0) . snd) obj))
+
+writeGLPLP :: (Show v, Ord v, Real c) => FilePath -> LP v c -> GLPK ()
+writeGLPLP file lp = do
+	vars <- writeProblem lp
+	sequence_ [setColName i (show v) | (v, i) <- assocs vars]
+	sequence_ [setRowName i lab | (i, Constr (Just lab) _ _) <- zip [1..] (constraints lp)]
+	writeLP file
diff --git a/src/Data/LinearProgram/GLPK/Internal.hs b/src/Data/LinearProgram/GLPK/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK/Internal.hs
@@ -0,0 +1,173 @@
+{-# LANGUAGE RecordWildCards, ScopedTypeVariables, ForeignFunctionInterface, BangPatterns #-}
+module Data.LinearProgram.GLPK.Internal (writeProblem, solveSimplex, mipSolve,
+	getObjVal, getRowPrim, getColPrim, mipObjVal, mipRowVal, mipColVal, getBadRay) where
+{-(writeProblem, addCols,
+	addRows, createIndex, findCol, findRow, getColPrim, getRowPrim, getObjVal,
+	mipColVal, mipRowVal, mipObjVal, mipSolve, setColBounds, setColKind, setColName, setMatRow,
+	setObjCoef, setObjectiveDirection, setRowBounds, setRowName, solveSimplex) where-}
+
+import Control.Monad
+import Prelude hiding ((+),(*))
+import Foreign.Ptr
+import Foreign.C
+import Foreign.Marshal.Array
+
+import Data.Bits
+import Data.Map hiding (map)
+-- import Data.Bounds
+import Data.LinearProgram.Common
+import Data.LinearProgram.GLPK.Types
+
+-- foreign import ccall "c_glp_set_obj_name" glpSetObjName :: Ptr GlpProb -> CString -> IO ()
+-- foreign import ccall unsafe "c_glp_set_obj_dir" glpSetObjDir :: Ptr GlpProb -> CInt -> IO ()
+foreign import ccall unsafe "c_glp_minimize" glpMinimize :: Ptr GlpProb -> IO ()
+foreign import ccall unsafe "c_glp_maximize" glpMaximize :: Ptr GlpProb -> IO ()
+foreign import ccall unsafe "c_glp_add_rows" glpAddRows :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_add_cols" glpAddCols :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_set_row_bnds" glpSetRowBnds :: Ptr GlpProb -> CInt -> CInt -> CDouble -> CDouble -> IO ()
+foreign import ccall unsafe "c_glp_set_col_bnds" glpSetColBnds :: Ptr GlpProb -> CInt -> CInt -> CDouble -> CDouble -> IO ()
+foreign import ccall unsafe "c_glp_set_obj_coef" glpSetObjCoef :: Ptr GlpProb -> CInt -> CDouble -> IO ()
+foreign import ccall unsafe "c_glp_set_mat_row" glpSetMatRow :: Ptr GlpProb -> CInt -> CInt -> Ptr CInt -> Ptr CDouble -> IO ()
+-- foreign import ccall unsafe "c_glp_create_index" glpCreateIndex :: Ptr GlpProb -> IO ()
+-- foreign import ccall unsafe "c_glp_find_row" glpFindRow :: Ptr GlpProb -> CString -> IO CInt
+-- foreign import ccall unsafe "c_glp_find_col" glpFindCol :: Ptr GlpProb -> CString -> IO CInt
+foreign import ccall unsafe "c_glp_solve_simplex" glpSolveSimplex :: Ptr GlpProb -> CInt -> CInt -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_get_obj_val" glpGetObjVal :: Ptr GlpProb -> IO CDouble
+foreign import ccall unsafe "c_glp_get_row_prim" glpGetRowPrim :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_get_col_prim" glpGetColPrim :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_set_col_kind" glpSetColKind :: Ptr GlpProb -> CInt -> CInt -> IO ()
+foreign import ccall unsafe "c_glp_mip_solve" glpMipSolve ::
+	Ptr GlpProb -> CInt -> CInt -> CInt -> CInt -> CInt -> CInt -> CInt -> CDouble -> CInt -> IO CInt
+foreign import ccall unsafe "c_glp_mip_obj_val" glpMIPObjVal :: Ptr GlpProb -> IO CDouble
+foreign import ccall unsafe "c_glp_mip_row_val" glpMIPRowVal :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_mip_col_val" glpMIPColVal :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall unsafe "c_glp_set_row_name" glpSetRowName :: Ptr GlpProb -> CInt -> CString -> IO ()
+foreign import ccall unsafe "c_glp_get_bad_ray" glpGetBadRay :: Ptr GlpProb -> IO CInt
+
+setObjectiveDirection :: Direction -> GLPK ()
+setObjectiveDirection dir = GLP $ case dir of
+	Min	-> glpMinimize
+	Max	-> glpMaximize
+
+getBadRay :: GLPK (Maybe Int)
+getBadRay = liftM (\ x -> guard (x /= 0) >> return (fromIntegral x)) $ GLP glpGetBadRay
+
+addRows :: Int -> GLPK Int
+addRows n = GLP $ liftM fromIntegral . flip glpAddRows (fromIntegral n)
+
+addCols :: Int -> GLPK Int
+addCols n = GLP $ liftM fromIntegral . flip glpAddCols (fromIntegral n)
+
+setRowBounds :: Real a => Int -> Bounds a -> GLPK ()
+setRowBounds i bds = GLP $ \ lp -> onBounds (glpSetRowBnds lp (fromIntegral i)) bds
+
+setColBounds :: Real a => Int -> Bounds a -> GLPK ()
+setColBounds i bds = GLP $ \ lp -> onBounds (glpSetColBnds lp (fromIntegral i)) bds
+
+onBounds :: Real a => (CInt -> CDouble -> CDouble -> x) -> Bounds a -> x
+onBounds f bds = case bds of
+	Free		-> f 1 0 0
+	LBound a	-> f 2 (realToFrac a) 0
+	UBound a	-> f 3 0 (realToFrac a)
+	Bound a b	-> f 4 (realToFrac a) (realToFrac b)
+	Equ a		-> f 5 (realToFrac a) 0
+
+{-# SPECIALIZE setObjCoef :: Int -> Double -> GLPK (), Int -> Int -> GLPK () #-}
+setObjCoef :: Real a => Int -> a -> GLPK ()
+setObjCoef i v = GLP $ \ lp -> glpSetObjCoef lp (fromIntegral i) (realToFrac v)
+
+{-# SPECIALIZE setMatRow :: Int -> [(Int, Double)] -> GLPK (), Int -> [(Int, Int)] -> GLPK () #-}
+setMatRow :: Real a => Int -> [(Int, a)] -> GLPK ()
+setMatRow i row = GLP $ \ lp ->
+	allocaArray (len+1) $ \ (ixs :: Ptr CInt) -> allocaArray (len+1) $ \ (coeffs :: Ptr CDouble) -> do
+		pokeArray ixs (0:map (fromIntegral . fst) row)
+		pokeArray coeffs (0:map (realToFrac . snd) row)
+		glpSetMatRow lp (fromIntegral i) (fromIntegral len) ixs coeffs
+	where	len = length row
+
+-- createIndex :: GLPK ()
+-- createIndex = GLP glpCreateIndex
+
+-- findRow :: String -> GLPK Int
+-- findRow nam = GLP $ liftM fromIntegral . withCString nam . glpFindRow
+
+-- findCol :: String -> GLPK Int
+-- findCol nam = GLP $ liftM fromIntegral . withCString nam . glpFindCol
+
+solveSimplex :: MsgLev -> Int -> Bool -> GLPK ReturnCode
+solveSimplex msglev tmLim presolve = GLP $ \ lp -> liftM (toEnum . fromIntegral) $ glpSolveSimplex lp
+	(getMsgLev msglev)
+	tmLim'
+	(if presolve then 1 else 0)
+	where	tmLim' = fromIntegral (tmLim * 1000)
+
+getMsgLev :: MsgLev -> CInt
+getMsgLev = fromIntegral . fromEnum
+
+getObjVal :: GLPK Double
+getObjVal = liftM realToFrac $ GLP glpGetObjVal
+
+getRowPrim :: Int -> GLPK Double
+getRowPrim i = liftM realToFrac $ GLP (`glpGetRowPrim` fromIntegral i)
+
+getColPrim :: Int -> GLPK Double
+getColPrim i = liftM realToFrac $ GLP (`glpGetColPrim` fromIntegral i)
+
+setColKind :: Int -> VarKind -> GLPK ()
+setColKind i kind = GLP $ \ lp -> glpSetColKind lp (fromIntegral i) (fromIntegral $ 1 + fromEnum kind)
+
+mipSolve :: MsgLev -> BranchingTechnique -> BacktrackTechnique -> Preprocessing -> Bool ->
+	[Cuts] -> Double -> Int -> Bool -> GLPK ReturnCode
+mipSolve msglev brt btt pp fp cuts mipgap tmlim presol =
+		liftM (toEnum . fromIntegral) $ GLP $ \ lp -> glpMipSolve lp msglev'
+						brt' btt' pp' fp' tmlim' cuts' mipgap' presol'
+	where	!msglev' = getMsgLev msglev
+		!brt' = 1 + fromIntegral (fromEnum brt)
+		!btt' = 1 + fromIntegral (fromEnum btt)
+		!pp' = fromIntegral (fromEnum pp)
+		!fp' = fromIntegral (fromEnum fp)
+		!cuts' = (if GMI `elem` cuts then 1 else 0) .|.
+			(if MIR `elem` cuts then 2 else 0) .|.
+			(if Cov `elem` cuts then 4 else 0) .|.
+			(if Clq `elem` cuts then 8 else 0)
+		!mipgap' = realToFrac mipgap
+		!tmlim' = fromIntegral (1000 * tmlim)
+		!presol' = fromIntegral (fromEnum presol)
+
+mipObjVal :: GLPK Double
+mipObjVal = liftM realToFrac $ GLP glpMIPObjVal
+
+mipRowVal :: Int -> GLPK Double
+mipRowVal i = liftM realToFrac $ GLP (`glpMIPRowVal` fromIntegral i)
+
+mipColVal :: Int -> GLPK Double
+mipColVal i = liftM realToFrac $ GLP (`glpMIPColVal` fromIntegral i)
+
+setRowName :: Int -> String -> GLPK ()
+setRowName i nam = GLP $ withCString nam . flip glpSetRowName (fromIntegral i)
+
+{-# SPECIALIZE writeProblem :: Ord v => LP v Double -> GLPK (Map v Int),
+	Ord v => LP v Int -> GLPK (Map v Int) #-}
+writeProblem :: (Ord v, Real c) => LP v c -> GLPK (Map v Int)
+writeProblem LP{..} = do
+	setObjectiveDirection direction
+	i0 <- addCols nVars
+	let allVars' = fmap (i0 +) allVars
+	sequence_ [setObjCoef i v | (i, v) <- elems $ intersectionWith (,) allVars' objective]
+	j0 <- addRows (length constraints)
+	sequence_ [do	maybe (return ()) (setRowName j) lab
+			setMatRow j
+				[(i, v) | (i, v) <- elems (intersectionWith (,) allVars' f)]
+			setRowBounds j bnds
+				| (j, Constr lab f bnds) <- zip [j0..] constraints]
+-- 	createIndex
+	sequence_ [setColBounds i bnds |
+			(i, bnds) <- elems $ intersectionWith (,) allVars' varBounds]
+	sequence_ [setColBounds i Free | i <- elems $ difference allVars' varBounds]
+	sequence_ [setColKind i knd |
+			(i, knd) <- elems $ intersectionWith (,) allVars' varTypes]
+	return allVars'
+	where	allVars0 = fmap (const ()) objective `union`
+			unions [fmap (const ()) f | Constr _ f _ <- constraints] `union`
+			fmap (const ()) varBounds `union` fmap (const ()) varTypes
+		(nVars, allVars) = mapAccum (\ n _ -> (n+1, n)) (0 :: Int) allVars0
diff --git a/src/Data/LinearProgram/GLPK/Solver.hs b/src/Data/LinearProgram/GLPK/Solver.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK/Solver.hs
@@ -0,0 +1,119 @@
+{-# OPTIONS -funbox-strict-fields #-}
+{-# LANGUAGE TupleSections, RecordWildCards #-}
+
+-- | Interface between the Haskell representation of a linear programming problem, a value of type 'LP', and
+-- the GLPK solver.  The options available to the solver correspond naturally with GLPK's available options,
+-- so to find the meaning of any particular option, consult the GLPK documentation.
+-- 
+-- The option of which solver to use -- the general LP solver, which solves a problem over the reals, or the 
+-- MIP solver, which allows variables to be restricted to integers -- can be made by choosing the appropriate
+-- constructor for 'GLPOpts'.
+-- 
+-- The marshalling from Haskell to C is specialized for 'Int's and 'Double's, so using those types in your
+-- linear program is recommended.
+module Data.LinearProgram.GLPK.Solver (
+	-- * Solver options
+	GLPOpts(..),
+	simplexDefaults, 
+	mipDefaults, 
+	-- * Running the solver
+	glpSolveVars,
+	RowValue(..),
+	glpSolveAll,
+	-- * GLPK enumerations
+	ReturnCode(..),
+	MsgLev(..), 
+	BranchingTechnique(..),
+	BacktrackTechnique(..), 
+	Preprocessing(..), 
+	Cuts(..)) where 
+
+import Control.Monad
+
+import Data.Map
+import Data.LinearProgram.Spec
+import Data.LinearProgram.GLPK.Common
+
+-- | Options available for customizing GLPK operations.  This also determines
+-- which kind of solving is performed -- relaxed LP, or MIP.
+data GLPOpts = SimplexOpts {msgLev :: MsgLev, tmLim :: !Int, presolve :: Bool} |
+	MipOpts {msgLev :: MsgLev, tmLim :: !Int, presolve :: Bool,
+		brTech :: BranchingTechnique, btTech :: BacktrackTechnique,
+		ppTech :: Preprocessing,
+		fpHeur :: Bool,
+		cuts :: [Cuts],
+		mipGap :: !Double}
+
+data RowValue v c = RowVal {row :: !(Constraint v c), rowVal :: !Double}
+
+simplexDefaults, mipDefaults :: GLPOpts
+simplexDefaults = SimplexOpts MsgOn 10000 True
+mipDefaults = MipOpts MsgOn 10000 True DrTom LocBound AllPre False [] 0.0
+
+{-# SPECIALIZE glpSolveVars :: Ord v => GLPOpts -> LP v Double -> IO (ReturnCode, Maybe (Double, Map v Double)),
+	Ord v => GLPOpts -> LP v Int -> IO (ReturnCode, Maybe (Double, Map v Double)) #-}
+-- | Solves the linear or mixed integer programming problem.  Returns
+-- the value of the objective function, and the values of the variables.
+glpSolveVars :: (Ord v, Real c) => GLPOpts -> LP v c -> IO (ReturnCode, Maybe (Double, Map v Double))
+glpSolveVars opts@SimplexOpts{} lp = runGLPK $ do
+	(code, vars) <- doGLP opts lp
+	liftM (code, ) $ maybe (return Nothing) ( \ vars -> do
+		obj <- getObjVal
+		vals <- sequence [do
+			val <- getColPrim i
+			return (v, val)
+				| (v, i) <- assocs vars]
+		return (Just (obj, fromDistinctAscList vals))) vars
+glpSolveVars opts@MipOpts{} lp = runGLPK $ do
+	(code, vars) <- doGLP opts lp
+	liftM (code, ) $ maybe (return Nothing) (\ vars -> do
+		obj <- mipObjVal
+		vals <- sequence [do
+			val <- mipColVal i
+			return (v, val)
+				| (v, i) <- assocs vars]
+		return (Just (obj, fromDistinctAscList vals))) vars
+
+{-# SPECIALIZE glpSolveAll :: 
+	Ord v => GLPOpts -> LP v Double -> IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v Double])),
+	Ord v => GLPOpts -> LP v Int -> IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v Int])) #-}
+-- | Solves the linear or mixed integer programming problem.  Returns
+-- the value of the objective function, the values of the variables,
+-- and the values of any labeled rows.
+glpSolveAll :: (Ord v, Real c) => GLPOpts -> LP v c -> IO (ReturnCode, Maybe (Double, Map v Double, [RowValue v c]))
+glpSolveAll opts@SimplexOpts{} lp@LP{..} = runGLPK $ do
+	(code, vars) <- doGLP opts lp
+	liftM (code, ) $ maybe (return Nothing) (\ vars -> do
+		obj <- getObjVal
+		vals <- sequence [do
+			val <- getColPrim i
+			return (v, val)
+				| (v, i) <- assocs vars]
+		rows <- sequence [liftM (RowVal c) (getRowPrim i)
+					| (i, c) <- zip [1..] constraints]
+		return (Just (obj, fromDistinctAscList vals, rows))) vars
+glpSolveAll opts@MipOpts{} lp@LP{..} = runGLPK $ do
+	(code, vars) <- doGLP opts lp
+	liftM (code, ) $ maybe (return Nothing) (\ vars -> do
+		obj <- mipObjVal
+		vals <- sequence [do
+			val <- mipColVal i
+			return (v, val)
+				| (v, i) <- assocs vars]
+		rows <- sequence [liftM (RowVal c) (mipRowVal i)
+					| (i, c) <- zip [1..] constraints]
+		return (Just (obj, fromDistinctAscList vals, rows))) vars
+
+{-# SPECIALIZE doGLP :: Ord v => GLPOpts -> LP v Double -> GLPK (ReturnCode, Maybe (Map v Int)),
+	Ord v => GLPOpts -> LP v Int -> GLPK (ReturnCode, Maybe (Map v Int)) #-}
+doGLP :: (Ord v, Real c) => GLPOpts -> LP v c -> GLPK (ReturnCode, Maybe (Map v Int))
+doGLP SimplexOpts{..} lp = do
+	vars <- writeProblem lp
+	success <- solveSimplex msgLev tmLim presolve
+	bad <- getBadRay
+	maybe (return (success, guard (gaveAnswer success) >> Just vars)) (fail . show) bad
+doGLP MipOpts{..} lp = do
+	vars <- writeProblem lp
+	success <- mipSolve msgLev brTech btTech ppTech fpHeur cuts mipGap tmLim presolve
+	bad <- getBadRay
+	return (success, guard (gaveAnswer success) >> Just vars)
diff --git a/src/Data/LinearProgram/GLPK/Types.hs b/src/Data/LinearProgram/GLPK/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/GLPK/Types.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE EmptyDataDecls, ForeignFunctionInterface #-}
+
+module Data.LinearProgram.GLPK.Types where
+
+import Control.Concurrent (runInBoundThread)
+import Control.Exception (bracket)
+import qualified Control.Monad.Fail as Fail
+import Control.Monad.Trans (MonadIO (..))
+import Control.Monad (ap)
+
+import Foreign.Ptr
+import Foreign.ForeignPtr
+
+foreign import ccall unsafe "c_glp_create_prob" glpCreateProb :: IO (Ptr GlpProb)
+foreign import ccall unsafe "c_glp_delete_prob" glpDelProb :: Ptr GlpProb -> IO ()
+
+data GlpProb
+
+data ReturnCode = Success | InvalidBasis | SingularMatrix | IllConditionedMatrix | 
+        InvalidBounds | SolverFailed | ObjLowerLimReached | ObjUpperLimReached | 
+        IterLimReached | TimeLimReached | NoPrimalFeasible | NoDualFeasible | RootLPOptMissing |
+        SearchTerminated | MipGapTolReached | NoPrimDualFeasSolution | NoConvergence |
+        NumericalInstability | InvalidData | ResultOutOfRange deriving (Eq, Show, Enum)
+
+gaveAnswer :: ReturnCode -> Bool
+gaveAnswer = flip elem [Success, IterLimReached, TimeLimReached, SearchTerminated, MipGapTolReached]
+
+newtype GLPK a = GLP {execGLPK :: Ptr GlpProb -> IO a}
+
+runGLPK :: GLPK a -> IO a
+runGLPK m = runInBoundThread $ bracket glpCreateProb glpDelProb (execGLPK m)
+
+instance Monad GLPK where
+        {-# INLINE return #-}
+        {-# INLINE (>>=) #-}
+        return x = GLP $ \ _ -> return x
+        m >>= k = GLP $ \ lp -> do      x <- execGLPK m lp
+                                        execGLPK (k x) lp
+instance Fail.MonadFail GLPK where
+        fail s = GLP . const $ Fail.fail s
+instance Functor GLPK where
+  fmap f (GLP k) = GLP $ \p -> fmap f (k p)
+
+instance Applicative GLPK where
+  pure = return
+  (<*>) = ap
+
+instance MonadIO GLPK where
+        liftIO m = GLP (const m)
+
+data MsgLev = MsgOff | MsgErr | MsgOn | MsgAll deriving (Eq, Enum, Read, Show)
+data BranchingTechnique = FirstFrac | LastFrac | MostFrac | DrTom | HybridP deriving (Eq, Enum, Read, Show)
+data BacktrackTechnique = DepthFirst | BreadthFirst | LocBound | ProjHeur deriving (Eq, Enum, Read, Show)
+data Preprocessing = NoPre | RootPre | AllPre deriving (Eq, Enum, Read, Show)
+data Cuts = GMI | MIR | Cov | Clq deriving (Eq, Enum, Read, Show)
diff --git a/src/Data/LinearProgram/LinExpr.hs b/src/Data/LinearProgram/LinExpr.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/LinExpr.hs
@@ -0,0 +1,62 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+module Data.LinearProgram.LinExpr (LinExpr(..), solve, substituteExpr, simplifyExpr,
+	constTerm, coeffTerm, funcToExpr) where
+import Control.Monad
+
+import Data.LinearProgram.Types
+import Algebra.Classes
+import Algebra.Linear
+import Data.Functor
+import Data.Foldable
+
+import Data.Map
+
+import Prelude hiding (lookup, filter, foldr, Num(..), recip)
+
+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, Additive c) => Additive (LinExpr v c) where
+	zero = LinExpr zero zero
+	LinExpr a1 c1 + LinExpr a2 c2 = LinExpr (a1 + a2) (c1 + c2)
+
+instance (Ord v, Group c) => Group (LinExpr v c) where
+	LinExpr a1 c1 - LinExpr a2 c2 = LinExpr (a1 - a2) (c1 - c2)
+	negate (LinExpr a c) = LinExpr (negate a) (negate c)
+
+instance (Ord v,AbelianAdditive c) => AbelianAdditive (LinExpr v c)
+
+instance (Ord v, Ring c) => Module c (LinExpr v c) where
+	k *^ LinExpr a c = LinExpr (k *^ a) (k * c)
+
+substituteExpr :: (Ord v, Ring 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, Ring 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 (negate 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 = negate (recip 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/src/Data/LinearProgram/Spec.hs b/src/Data/LinearProgram/Spec.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/Spec.hs
@@ -0,0 +1,155 @@
+{-# LANGUAGE TupleSections, RecordWildCards, DeriveFunctor #-}
+module Data.LinearProgram.Spec (Constraint(..), VarTypes, ObjectiveFunc, VarBounds, LP(..),
+        mapVars, mapVals, allVars, linCombination) where
+
+import Prelude hiding (negate, (+))
+import Control.DeepSeq
+import Control.Monad
+import Data.Char (isSpace)
+import Data.Map hiding (map, foldl)
+
+import Text.ParserCombinators.ReadP
+
+import Algebra.Classes
+import Data.LinearProgram.Types
+import qualified Data.Map as M
+
+-- | Representation of a linear constraint on the variables, possibly labeled.
+-- The function may be bounded both above and below.
+data Constraint v c = Constr (Maybe String)
+                        (LinFunc v c)
+                        (Bounds c) deriving (Functor)
+-- | A mapping from variables to their types.  Variables not mentioned are assumed to be continuous,
+type VarTypes v = Map v VarKind
+-- | An objective function for a linear program.
+type ObjectiveFunc = LinFunc
+-- | A mapping from variables to their boundaries.  Variables not mentioned are assumed to be free.
+type VarBounds v c = Map v (Bounds c)
+
+-- | The specification of a linear programming problem with variables in @v@ and coefficients/constants in @c@.
+--   Note: the 'Read' and 'Show' implementations do not correspond to any particular linear program specification format.
+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)
+
+linCombination :: (Ord v, Additive r) => [(r, v)] -> LinFunc v r
+linCombination xs = M.fromListWith (+) [(v, r) | (r, v) <- xs]
+
+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"
+        Equ x   -> expr ++ " = " ++ show x
+        LBound x -> expr ++ " >= " ++ show x
+        UBound x -> expr ++ " <= " ++ show x
+        Bound l u -> show l ++ " <= " ++ expr ++ " <= " ++ show u
+
+showFunc :: (Show v, Ord c, Show c, Num c, Group c) => LinFunc v c -> String
+showFunc func = case assocs func of
+        []      -> "0"
+        ((v,c):vcs) ->
+                show c ++ " " ++ map replaceSpace (show v) ++
+                        concatMap showTerm vcs
+        where   showTerm (v, c) = case compare c 0 of
+                        EQ      -> ""
+                        GT      -> " + " ++ show c ++ " " ++ show v
+                        LT      -> " - " ++ show (negate c) ++ " " ++ show v
+
+replaceSpace :: Char -> Char
+replaceSpace c
+        | isSpace c     = '_'
+        | otherwise     = c
+
+instance (Show v, Ord c, Show c, Num c, Group c) => Show (Constraint v c) where
+        show (Constr lab func bds) = maybe "" (++ ": ") lab ++
+                showBds (showFunc func) bds
+
+instance (Read v, Ord v, Read c, Ord c, Num c, Group c) => Read (Constraint v c) where
+        readsPrec _= readP_to_S $ liftM toConstr (lab <++ nolab) where
+                toConstr (l, f, bds) = Constr l (fromList f) bds
+                lab = do        skipSpaces
+                                label <- manyTill get (skipSpaces >> char ':')
+                                (_, f, bds) <- nolab
+                                return (Just label, f, bds)
+                nolab = liftM (\ (f, bds) -> (Nothing, f, bds)) $ readBds readConst readFunc
+                readFunc = (do  c <- readCoef readConst
+                                v <- readVar
+                                liftM ((v, c):) readFunc) <++ return []
+                readConst = readS_to_P reads
+                readVar = readS_to_P reads
+
+readCoef :: (Num c, Group c) => ReadP c -> ReadP c
+readCoef readC = between skipSpaces skipSpaces $
+        (do     char '+'
+                skipSpaces
+                readC') <++
+        (do     char '-'
+                skipSpaces
+                negate <$> readC') <++ readC'
+        where   readC' = readC <++ return 1
+
+optMaybe :: ReadP a -> ReadP (Maybe a)
+optMaybe p = fmap Just p <++ return Nothing
+
+readBds :: Ord c => ReadP c -> ReadP a -> ReadP (a, Bounds c)
+readBds cst expr = do
+        left <- optMaybe (do    lb <- cst
+                                skipSpaces
+                                rel <- readRelation
+                                return (lb, rel))
+        skipSpaces
+        f <- expr
+        skipSpaces
+        right <- optMaybe (do   rel <- readRelation
+                                skipSpaces
+                                ub <- cst
+                                return (ub, revOrd rel))
+        return (f, getBd left `mappend` getBd right)
+        where   revOrd :: Ordering -> Ordering
+                revOrd GT = LT
+                revOrd LT = GT
+                revOrd EQ = EQ
+                getBd :: Maybe (c, Ordering) -> Bounds c
+                getBd Nothing = Free
+                getBd (Just (x, cmp)) = case cmp of
+                        EQ      -> Equ x
+                        GT      -> LBound x
+                        LT      -> UBound x
+                readRelation = choice [char '<' >> optional (char '=') >> return LT,
+                        char '=' >> return EQ,
+                        char '>' >> optional (char '=') >> return GT]
+
+{-# SPECIALIZE mapVars :: Ord v' => (v -> v') -> LP v Double -> LP v' Double #-}
+-- | Applies the specified function to the variables in the linear program.
+-- If multiple variables in the original program are mapped to the same variable in the new program,
+-- in general, we set those variables to all be equal, as follows.
+--
+-- * In linear functions, including the objective function and the constraints,
+--      coefficients will be added together.  For instance, if @v1,v2@ are mapped to the same
+--      variable @v'@, then a linear function of the form @c1 *& v1 ^+^ c2 *& v2@ will be mapped to
+--      @(c1 ^+^ c2) *& v'@.
+--
+-- * In variable bounds, bounds will be combined.  An error will be thrown if the bounds
+--      are mutually contradictory.
+--
+-- * In variable kinds, the most restrictive kind will be retained.
+mapVars :: (Ord v', Ord c, Group c) => (v -> v') -> LP v c -> LP v' c
+mapVars f LP{..} =
+        LP{objective = mapKeysWith (+) f objective,
+                constraints = [Constr lab (mapKeysWith (+) f func) bd | Constr lab func bd <- constraints],
+                varBounds = mapKeysWith mappend f varBounds,
+                varTypes = mapKeysWith mappend f varTypes, ..}
+
+-- | Applies the specified function to the constants in the linear program.  This is only safe
+-- for a monotonic function.
+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 (LP v c) where
+        rnf LP{..} = direction `deepseq` objective `deepseq` constraints `deepseq`
+                varBounds `deepseq` rnf varTypes
diff --git a/src/Data/LinearProgram/Types.hs b/src/Data/LinearProgram/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearProgram/Types.hs
@@ -0,0 +1,82 @@
+{-# LANGUAGE DeriveFunctor, DeriveGeneric #-}
+module Data.LinearProgram.Types (LinFunc, VarKind(..), Direction(..), Bounds(..)) where
+
+import Control.DeepSeq
+import Data.Monoid
+import GHC.Generics
+import Data.Map
+
+type LinFunc = Map
+
+
+data VarKind = ContVar | IntVar | BinVar deriving (Eq, Ord, Enum, Show, Read, Generic)
+
+-- instance NFData VarKind
+
+instance Semigroup VarKind where
+        (<>) = max
+
+instance Monoid VarKind where
+        mempty = ContVar
+        mappend = (<>)
+
+data Direction = Min | Max deriving (Eq, Ord, Enum, Show, Read, Generic)
+
+-- instance NFData Direction
+
+data Bounds a =
+        Free | LBound !a | UBound !a | Equ !a | Bound !a !a deriving (Eq, Show, Read, Functor)
+
+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 (Bounds a)
+
+-- Bounds form a monoid under intersection.
+instance Ord a => Monoid (Bounds a) where
+        mempty = Free
+        mappend = (<>)
+
+instance Ord a => Semigroup (Bounds a) where
+        Free <> bd = bd
+        bd <> Free = bd
+        Equ a <> Equ b
+                | a == b        = Equ a
+        Equ a <> UBound b
+                | a <= b        = Equ a
+        Equ a <> LBound b
+                | a >= b        = Equ a
+        Equ a <> Bound l u
+                | a >= l && a <= u
+                                = Equ a
+        Equ _ <> _ = infeasible
+        UBound b <> Equ a
+                | a <= b        = Equ a
+        LBound b <> Equ a
+                | a >= b        = Equ a
+        Bound l u <> Equ a
+                | a >= l && a <= u
+                                = Equ a
+        _ <> Equ _ = infeasible
+        LBound a <> LBound b = LBound (max a b)
+        LBound l <> UBound u = bound l u
+        UBound u <> LBound l = bound l u
+        LBound a <> Bound l u = bound (max a l) u
+        Bound l u <> LBound a = bound (max a l) u
+        UBound a <> UBound b = UBound (min a b)
+        UBound a <> Bound l u = bound l (min a u)
+        Bound l u <> UBound a = bound l (min a u)
+        Bound l u <> Bound l' u' = bound (max l l') (min u u')
+
+infeasible :: Bounds a
+infeasible = error "Mutually contradictory constraints found."
+
+bound :: Ord a => a -> a -> Bounds a
+bound l u       | l <= u        = Bound l u
+                | otherwise     = infeasible