diff --git a/Data/LinFunc.hs b/Data/LinFunc.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinFunc.hs
@@ -0,0 +1,121 @@
+{-# LANGUAGE UndecidableInstances, FlexibleInstances, MultiParamTypeClasses #-}
+
+module Data.LinFunc (LinFunc, Module(..), var, varSum, (*&), vsum, combination, linCombination) where
+
+import Control.Monad
+
+import qualified Data.Map as M
+import qualified Data.IntMap as IM
+import Data.Ratio
+import Data.Array.Base
+import Data.Array.IArray
+import Data.Array.Unboxed
+
+-- import Data.LinFunc.TH
+import Data.LinFunc.Class
+
+-- | @'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
+
+instance Module Int Int where
+	(*^) = (*)
+	zero = 0
+	(^+^) = (+)
+	(^-^) = (-)
+	neg = negate
+
+instance Module Double Double where
+	(*^) = (*)
+	zero = 0
+	(^+^) = (+)
+	(^-^) = (-)
+	neg = negate
+
+instance Module Integer Integer where
+	(*^) = (*)
+	zero = 0
+	(^+^) = (+)
+	(^-^) = (-)
+	neg = negate
+
+
+instance Integral a => Module (Ratio a) (Ratio a) where
+	{-# SPECIALIZE instance Module Rational Rational #-}
+	{-# SPECIALIZE instance Module (Ratio Int) (Ratio Int) #-}
+	(*^) = (*)
+	zero = 0
+	(^+^) = (+)
+	(^-^) = (-)
+	neg = negate
+
+instance Module r m => Module r (a -> m) where
+	(*^) = fmap . (*^)
+	zero = const zero
+	(^+^) = liftM2 (^+^)
+	(^-^) = liftM2 (^-^)
+	neg = fmap neg
+
+instance (Ord k, Module r m) => Module r (M.Map k m) where
+	(*^) = fmap . (*^)
+	zero = M.empty
+	(^+^) = M.unionWith (^+^)
+	neg = fmap neg
+
+instance Module r m => Module r (IM.IntMap m) where
+	(*^) = fmap . (*^)
+	zero = IM.empty
+	(^+^) = IM.unionWith (^+^)
+	neg = fmap neg
+	
+instance (Module r m) => Module r (Array Int m) where
+	(*^) = amap . (*^)
+	zero = listArray (0,0) [zero]
+	a ^+^ b	| numElements a >= numElements b
+			= accum (^+^) a (assocs b)
+		| otherwise
+			= accum (^+^) b (assocs a)
+	a ^-^ b | numElements a >= numElements b
+			= accum (^-^) a (assocs b)
+		| otherwise
+			= accum (^-^) b (assocs a)
+	neg = amap neg
+
+instance (IArray UArray m, Module r m) => Module r (UArray Int m) where
+	(*^) = amap . (*^)
+	zero = listArray (0,0) [zero]
+	a ^+^ b	| numElements a >= numElements b
+			= accum (^+^) a (assocs b)
+		| otherwise
+			= accum (^+^) b (assocs a)
+	a ^-^ b | numElements a >= numElements b
+			= accum (^-^) a (assocs b)
+		| otherwise
+			= accum (^-^) b (assocs a)
+	neg = amap neg
+
+-- | Given a variable @v@, returns the function equivalent to @v@.
+var :: (Ord v, Num c) => v -> LinFunc v c
+var v = M.singleton v 1
+
+-- | @c '*&' v@ is equivalent to @c '*^' 'var' v@.
+(*&) :: (Ord v, Num c) => c -> v -> LinFunc v c
+c *& v = M.singleton v c
+
+-- | Equivalent to @'vsum' . 'map' 'var'@.
+varSum :: (Ord v, Num c) => [v] -> LinFunc v c
+varSum vs = M.fromList [(v, 1) | v <- vs]
+
+-- | Returns a vector sum.
+vsum :: Module r v => [v] -> v
+vsum = foldr (^+^) zero
+
+-- | Given a collection of vectors and scaling coefficients, returns this
+-- linear combination.
+combination :: Module r m => [(r, m)] -> m
+combination xs = vsum [r *^ m | (r, m) <- xs]
+
+-- | 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]
diff --git a/Data/LinFunc/Class.hs b/Data/LinFunc/Class.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinFunc/Class.hs
@@ -0,0 +1,19 @@
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}
+
+module Data.LinFunc.Class where
+
+infixr 4 ^+^
+infixr 4 ^-^
+infixr 6 *^
+
+-- | In algebra, if @r@ is a ring, an @r@-module is an additive group with a scalar multiplication
+-- operation.  When @r@ is a field, this is equivalent to a vector space.
+class Module r m | m -> r where
+	(*^) :: r -> m -> m
+	zero :: m
+	(^+^) :: m -> m -> m
+	(^-^) :: m -> m -> m
+	neg :: m -> m
+	
+	a ^-^ b = a ^+^ neg b
+	neg a = zero ^-^ a
diff --git a/Data/LinearProgram.hs b/Data/LinearProgram.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram.hs
@@ -0,0 +1,6 @@
+module Data.LinearProgram (module Data.LinearProgram.Spec, module Data.LinearProgram.Types,
+	module Data.LinFunc) where
+
+import Data.LinearProgram.Spec
+import Data.LinearProgram.Types
+import Data.LinFunc
diff --git a/Data/LinearProgram/GLPK.hs b/Data/LinearProgram/GLPK.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/GLPK.hs
@@ -0,0 +1,126 @@
+{-# OPTIONS -funbox-strict-fields #-}
+{-# LANGUAGE RecordWildCards #-}
+
+module Data.LinearProgram.GLPK (GLPOpts(..), MsgLev(..), BranchingTechnique(..),
+	BacktrackTechnique(..), Preprocessing(..), Cuts(..), 
+	simplexDefaults, mipDefaults, glpSolveVars, glpSolveAll) where
+
+import Control.Monad.Trans
+
+import Data.Map
+import Data.Maybe (catMaybes)
+import Data.LinearProgram.Spec
+import Data.LinearProgram.Types
+import Data.LinearProgram.GLPK.Internal
+
+import System.CPUTime
+
+import GHC.Exts(build)
+
+-- | 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}
+
+simplexDefaults, mipDefaults :: GLPOpts
+simplexDefaults = SimplexOpts MsgOn 10000 True
+mipDefaults = MipOpts MsgOn 10000 True DrTom LocBound AllPre False [] 0.0
+
+-- | 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 (Double, Map v Double)
+glpSolveVars opts@SimplexOpts{} lp = runGLPK $ do
+	Just vars <- doGLP opts lp
+	obj <- getObjVal
+	vals <- sequence [do
+		val <- getColPrim i
+		return (v, val)
+			| (v, i) <- assocs vars]
+	return (obj, fromDistinctAscList vals)
+glpSolveVars opts@MipOpts{} lp = runGLPK $ do
+	Just vars <- doGLP opts lp
+	obj <- mipObjVal
+	vals <- sequence [do
+		val <- mipColVal i
+		return (v, val)
+			| (v, i) <- assocs vars]
+	return (obj, fromDistinctAscList vals)
+
+-- | 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 (Double, Map v Double, Map String Double)
+glpSolveAll opts@SimplexOpts{} lp@LP{..} = runGLPK $ do
+	Just vars <- doGLP opts lp
+	obj <- getObjVal
+	vals <- sequence [do
+		val <- getColPrim i
+		return (v, val)
+			| (v, i) <- assocs vars]
+	rows <- sequence [maybe (return Nothing) (\ nam -> do
+				val <- getRowPrim i
+				return (Just (nam, val))) nam
+				| (i, Constr nam _ _) <- zip [0..] constraints]
+	return (obj, fromDistinctAscList vals, fromDistinctAscList (catMaybes rows))
+glpSolveAll opts@MipOpts{} lp@LP{..} = runGLPK $ do
+	Just vars <- doGLP opts lp
+	obj <- mipObjVal
+	vals <- sequence [do
+		val <- mipColVal i
+		return (v, val)
+			| (v, i) <- assocs vars]
+	rows <- sequence [maybe (return Nothing) (\ nam -> do
+				val <- mipRowVal i
+				return (Just (nam, val))) nam
+				| (i, Constr nam _ _) <- zip [0..] constraints]
+	return (obj, fromDistinctAscList vals, fromDistinctAscList (catMaybes rows))
+
+doGLP :: (Ord v, Real c) => GLPOpts -> LP v c -> GLPK (Maybe (Map v Int))
+doGLP SimplexOpts{..} lp = do
+	vars <- writeProblem lp
+	success <- solveSimplex msgLev tmLim presolve
+	return (if success then Just vars else Nothing)
+doGLP MipOpts{..} lp = do
+	vars <- writeProblem lp
+	time <- getTime
+	solveSimplex msgLev tmLim presolve
+	time' <- getTime
+	let tmLim' = (fromIntegral tmLim - time' + time + 1000000000000 - 1) `quot` 1000000000000
+	success <- mipSolve msgLev brTech btTech ppTech fpHeur cuts mipGap (fromIntegral tmLim') presolve
+	return (if success then Just vars else Nothing)
+	where	getTime = liftIO getCPUTime
+
+writeProblem :: (Ord v, Real c) => LP v c -> GLPK (Map v Int)
+writeProblem LP{..} = do
+	setObjectiveDirection direction
+	i0 <- addCols nVars
+	sequence_ [setObjCoef (i + i0) v | (i, v) <- elems $ intersectionWith (,) allVars objective]
+	j0 <- addRows (length constraints)
+	sequence_ [do	case lab of
+				Nothing	-> return ()
+				Just n	-> setRowName (j0 + j) n
+			setMatRow (j0 + j)
+				(elems (intersectionWith (,) allVars f))
+			setRowBounds (j0 + j) bnds
+				| (j, Constr lab f bnds) <- zip [0..] constraints]
+	createIndex
+	sequence_ [setColBounds (i0 + i) bnds |
+			(i, bnds) <- elems $ intersectionWith (,) allVars varBounds]
+	sequence_ [setColKind (i0 + 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
+		
+{-# 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/Internal.hs b/Data/LinearProgram/GLPK/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/GLPK/Internal.hs
@@ -0,0 +1,191 @@
+{-# LANGUAGE ScopedTypeVariables, EmptyDataDecls, ForeignFunctionInterface #-}
+module Data.LinearProgram.GLPK.Internal (GLPK, MsgLev (..), Preprocessing (..), Direction(..), BacktrackTechnique(..),
+	BranchingTechnique(..), Cuts(..), runGLPK, 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 Control.Monad.Trans
+
+import Debug.Trace
+
+import Foreign.Ptr
+import Foreign.C
+import Foreign.ForeignPtr
+import Foreign.Marshal.Array
+
+import Data.Bits
+-- import Data.Bounds
+import Data.LinearProgram.Types
+
+data GlpProb
+
+foreign import ccall "c_glp_create_prob" glpCreateProb :: IO (Ptr GlpProb)
+-- foreign import ccall "c_glp_set_obj_name" glpSetObjName :: Ptr GlpProb -> CString -> IO ()
+foreign import ccall "c_glp_set_obj_dir" glpSetObjDir :: Ptr GlpProb -> CInt -> IO ()
+foreign import ccall "c_glp_add_rows" glpAddRows :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall "c_glp_add_cols" glpAddCols :: Ptr GlpProb -> CInt -> IO CInt
+foreign import ccall "c_glp_set_row_name" glpSetRowName :: Ptr GlpProb -> CInt -> CString -> IO ()
+foreign import ccall "c_glp_set_col_name" glpSetColName :: Ptr GlpProb -> CInt -> CString -> IO ()
+foreign import ccall "c_glp_set_row_bnds" glpSetRowBnds :: Ptr GlpProb -> CInt -> CInt -> CDouble -> CDouble -> IO ()
+foreign import ccall "c_glp_set_col_bnds" glpSetColBnds :: Ptr GlpProb -> CInt -> CInt -> CDouble -> CDouble -> IO ()
+foreign import ccall "c_glp_set_obj_coef" glpSetObjCoef :: Ptr GlpProb -> CInt -> CDouble -> IO ()
+foreign import ccall "c_glp_set_mat_row" glpSetMatRow :: Ptr GlpProb -> CInt -> CInt -> Ptr CInt -> Ptr CDouble -> IO ()
+foreign import ccall "c_glp_delete_prob" glpDelProb :: Ptr GlpProb -> IO ()
+foreign import ccall "c_glp_create_index" glpCreateIndex :: Ptr GlpProb -> IO ()
+foreign import ccall "c_glp_find_row" glpFindRow :: Ptr GlpProb -> CString -> IO CInt
+foreign import ccall "c_glp_find_col" glpFindCol :: Ptr GlpProb -> CString -> IO CInt
+foreign import ccall "c_glp_solve_simplex" glpSolveSimplex :: Ptr GlpProb -> CInt -> CInt -> CInt -> IO CInt
+foreign import ccall "c_glp_get_obj_val" glpGetObjVal :: Ptr GlpProb -> IO CDouble
+foreign import ccall "c_glp_get_row_prim" glpGetRowPrim :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall "c_glp_get_col_prim" glpGetColPrim :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall "c_glp_set_col_kind" glpSetColKind :: Ptr GlpProb -> CInt -> CInt -> IO ()
+foreign import ccall "c_glp_mip_solve" glpMipSolve :: 
+	Ptr GlpProb -> CInt -> CInt -> CInt -> CInt -> CInt -> CInt -> CInt -> CDouble -> CInt -> IO CInt
+foreign import ccall "c_glp_mip_obj_val" glpMIPObjVal :: Ptr GlpProb -> IO CDouble
+foreign import ccall "c_glp_mip_row_val" glpMIPRowVal :: Ptr GlpProb -> CInt -> IO CDouble
+foreign import ccall "c_glp_mip_col_val" glpMIPColVal :: Ptr GlpProb -> CInt -> IO CDouble
+
+newtype GLPK a = GLP {execGLPK :: Ptr GlpProb -> IO a}
+
+runGLPK :: GLPK a -> IO a
+runGLPK m = do	lp <- glpCreateProb
+		ans <- execGLPK m lp
+		glpDelProb lp
+		return ans
+
+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 MonadIO GLPK where
+	liftIO m = GLP (const m)
+
+setObjectiveDirection :: Direction -> GLPK ()
+setObjectiveDirection dir = GLP $ flip glpSetObjDir 
+	(case dir of	Min	-> 1
+			Max	-> 2)
+
+addRows :: Int -> GLPK Int
+addRows n = GLP $ liftM (subtract 1 . fromIntegral) . flip glpAddRows (fromIntegral n)
+
+addCols :: Int -> GLPK Int
+addCols n = GLP $ liftM (subtract 1 . fromIntegral) . flip glpAddCols (fromIntegral n)
+
+setRowName :: Int -> String -> GLPK ()
+setRowName i nam = GLP $ withCString nam . flip glpSetRowName (fromIntegral (i+1))
+
+setColName :: Int -> String -> GLPK ()
+setColName i nam = GLP $ withCString nam . flip glpSetColName (fromIntegral (i+1))
+
+setRowBounds :: Real a => Int -> Bounds a -> GLPK ()
+setRowBounds i bds = GLP $ \ lp -> onBounds (glpSetRowBnds lp (fromIntegral (i+1))) bds
+
+setColBounds :: Real a => Int -> Bounds a -> GLPK ()
+setColBounds i bds = GLP $ \ lp -> onBounds (glpSetColBnds lp (fromIntegral (i+1))) 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
+
+setObjCoef :: Real a => Int -> a -> GLPK ()
+setObjCoef i v = GLP $ \ lp -> glpSetObjCoef lp (fromIntegral (i + 1)) (realToFrac v)
+
+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 . (+1) . fst) row)
+		pokeArray coeffs (0:map (realToFrac . snd) row)
+		glpSetMatRow lp (fromIntegral (i+1)) (fromIntegral len) ixs coeffs
+	where	len = length row
+
+createIndex :: GLPK ()
+createIndex = GLP glpCreateIndex
+
+findRow :: String -> GLPK Int
+findRow nam = GLP $ liftM (subtract 1 . fromIntegral) . withCString nam . glpFindRow
+
+findCol :: String -> GLPK Int
+findCol nam = GLP $ liftM (subtract 1 . fromIntegral) . withCString nam . glpFindCol
+
+data MsgLev = MsgOff | MsgErr | MsgOn | MsgAll
+
+solveSimplex :: MsgLev -> Int -> Bool -> GLPK Bool
+solveSimplex msglev tmLim presolve = GLP $ \ lp -> liftM (== 0) $ glpSolveSimplex lp
+	(getMsgLev msglev)
+	tmLim'
+	(if presolve then 1 else 0)
+	where	tmLim' = fromIntegral (tmLim * 1000)
+
+getMsgLev :: MsgLev -> CInt
+getMsgLev msglev = case msglev of
+	MsgOff	-> 0
+	MsgErr	-> 1
+	MsgOn	-> 2
+	MsgAll	-> 3
+
+getObjVal :: GLPK Double
+getObjVal = liftM realToFrac $ GLP glpGetObjVal
+
+getRowPrim :: Int -> GLPK Double
+getRowPrim i = liftM realToFrac $ GLP (`glpGetRowPrim` fromIntegral (i+1))
+
+getColPrim :: Int -> GLPK Double
+getColPrim i = liftM realToFrac $ GLP (`glpGetColPrim` fromIntegral (i+1))
+
+setColKind :: Int -> VarKind -> GLPK ()
+setColKind i kind = GLP $ \ lp -> glpSetColKind lp (fromIntegral (i+1)) (case kind of
+	ContVar -> 1
+	IntVar	-> 2
+	BinVar	-> 3)
+
+data BranchingTechnique = FirstFrac | LastFrac | MostFrac | DrTom | HybridP
+data BacktrackTechnique = DepthFirst | BreadthFirst | LocBound | ProjHeur
+data Preprocessing = NoPre | RootPre | AllPre
+data Cuts = GMI | MIR | Cov | Clq deriving (Eq)
+
+mipSolve :: MsgLev -> BranchingTechnique -> BacktrackTechnique -> Preprocessing -> Bool ->
+	[Cuts] -> Double -> Int -> Bool -> GLPK Bool
+mipSolve msglev brt btt pp fp cuts mipgap tmlim presol =
+		liftM (== 0) $ GLP $ \ lp -> glpMipSolve lp (getMsgLev msglev)
+						brt' btt' pp' fp' tmlim' cuts' mipgap' presol'
+	where	brt' = case brt of
+			FirstFrac	-> 1
+			LastFrac	-> 2
+			MostFrac	-> 3
+			DrTom		-> 4
+			HybridP		-> 5
+		btt' = case btt of
+			DepthFirst	-> 1
+			BreadthFirst	-> 2
+			LocBound	-> 3
+			ProjHeur	-> 4
+		pp' = case pp of
+			NoPre	-> 0
+			RootPre	-> 1
+			AllPre	-> 2
+		fp' = if fp then 1 else 0
+		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' = if presol then 1 else 0
+
+mipObjVal :: GLPK Double
+mipObjVal = liftM realToFrac $ GLP glpMIPObjVal
+
+mipRowVal :: Int -> GLPK Double
+mipRowVal i = liftM realToFrac $ GLP (`glpMIPRowVal` fromIntegral (i+1))
+
+mipColVal :: Int -> GLPK Double
+mipColVal i = liftM realToFrac $ GLP (`glpMIPRowVal` fromIntegral (i+1))
diff --git a/Data/LinearProgram/LPMonad.hs b/Data/LinearProgram/LPMonad.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/LPMonad.hs
@@ -0,0 +1,88 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module Data.LinearProgram.LPMonad where
+
+import Control.Monad.State.Strict
+
+import Data.Map
+import Data.Monoid
+-- import Data.Bounds
+
+import Data.LinFunc
+import Data.LinearProgram.Types
+import Data.LinearProgram.Spec
+
+-- | A 'State' monad used for the construction of a linear program.
+type LPM v c = State (LP v c)
+
+-- | Constructs a linear programming problem, returning any
+-- desired return value.
+runLPM :: (Ord v, Module r c) => LPM v c a -> (a, LP v c)
+runLPM m = runState m (LP Max zero [] mempty mempty)
+
+-- | Constructs a linear programming problem.
+execLPM :: (Ord v, Module r c) => LPM v c a -> LP v c
+execLPM = snd . runLPM
+
+-- | Sets the optimization direction of the linear program:
+-- maximization or minimization.
+setDirection :: Direction -> LPM v c ()
+setDirection dir = modify (\ lp -> lp{direction = dir})
+
+equal, leq, geq :: (Ord v, Module r c) => LinFunc v c -> LinFunc v c -> LPM v c ()
+equal f g = equalTo (f ^-^ g) zero
+leq f g = leqTo (f ^-^ g) zero
+geq = flip leq
+
+equal', leq', geq' :: (Ord v, Module r c) => String -> LinFunc v c -> LinFunc v c -> LPM v c ()
+equal' lab f g = equalTo' lab (f ^-^ g) zero
+leq' lab f g = leqTo' lab (f ^-^ g) zero
+geq' = flip . leq'
+
+equalTo, leqTo, geqTo :: LinFunc v c -> c -> LPM v c ()
+equalTo f v = constrain f (Equ v)
+leqTo f v = constrain f (UBound v)
+geqTo f v = constrain f (LBound v)
+
+equalTo', leqTo', geqTo' :: String -> LinFunc v c -> c -> LPM v c ()
+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)
+
+varEq, varLeq, varGeq :: (Ord v, Ord c) => v -> c -> LPM v c ()
+varEq v c = setVarBounds v (Equ c)
+varLeq v c = setVarBounds v (UBound c)
+varGeq v c = setVarBounds v (LBound c)
+
+varBds :: (Ord v, Ord c) => v -> c -> c -> LPM v c ()
+varBds v l u = setVarBounds v (Bound l u)
+
+constrain :: LinFunc v c -> Bounds c -> LPM v c ()
+constrain f bds = modify addConstr where
+	addConstr lp@LP{..}
+		= lp{constraints = Constr Nothing f bds:constraints}
+
+constrain' :: String -> LinFunc v c -> Bounds c -> LPM v c ()
+constrain' lab f bds = modify addConstr where
+	addConstr lp@LP{..}
+		= lp{constraints = Constr (Just lab) f bds:constraints}
+
+setObjective :: LinFunc v c -> LPM v c ()
+setObjective obj = modify setObj where
+	setObj lp = lp{objective = obj}
+
+addObjective :: (Ord v, Module r c) => LinFunc v c -> LPM v c ()
+addObjective obj = modify addObj where
+	addObj lp@LP{..}
+		= lp {objective = obj ^+^ objective}
+		
+addWeightedObjective :: (Ord v, Module r c) => r -> LinFunc v c -> LPM v c ()
+addWeightedObjective wt obj = addObjective (wt *^ obj)
+
+setVarBounds :: (Ord v, Ord c) => v -> Bounds c -> LPM v c ()
+setVarBounds var bds = modify addBds where
+	addBds lp@LP{..} = lp{varBounds = insertWith mappend var bds varBounds}
+
+setVarKind :: Ord v => v -> VarKind -> LPM v c ()
+setVarKind v k = modify setK where
+	setK lp@LP{..} = lp{varTypes = insertWith mappend v k varTypes}
diff --git a/Data/LinearProgram/Spec.hs b/Data/LinearProgram/Spec.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/Spec.hs
@@ -0,0 +1,16 @@
+module Data.LinearProgram.Spec where
+
+-- import Data.Bounds
+import Data.LinFunc
+import Data.LinearProgram.Types
+import Data.Map
+
+data Constraint v c = Constr (Maybe String)
+			(LinFunc v c)
+			(Bounds c) deriving (Read, Show)
+type VarTypes v = Map v VarKind
+type ObjectiveFunc = LinFunc
+type VarBounds v c = Map v (Bounds c)
+
+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)
diff --git a/Data/LinearProgram/Types.hs b/Data/LinearProgram/Types.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/Types.hs
@@ -0,0 +1,32 @@
+module Data.LinearProgram.Types where
+
+import Data.Monoid
+
+data VarKind = ContVar | IntVar | BinVar deriving (Eq, Ord, Show, Read)
+
+instance Monoid VarKind where
+	mempty = ContVar
+	mappend = max
+
+data Direction = Min | Max deriving (Eq, Ord, Show, Read)
+
+
+data Bounds a =
+	Free | LBound a | UBound a | Equ a | Bound a a deriving (Eq, Show, Read)
+
+-- 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 a
+	_ `mappend` Equ a = Equ a
+	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')
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,2 @@
+Copyright Louis Wasserman 2010
+GPL license
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/examples/example1.hs b/examples/example1.hs
new file mode 100644
--- /dev/null
+++ b/examples/example1.hs
@@ -0,0 +1,21 @@
+
+import Data.LinearProgram.LPMonad
+import Data.LinearProgram
+import Data.LinearProgram.GLPK
+
+objFun :: LinFunc String Int
+objFun = linCombination [(10, "x1"), (6, "x2"), (4, "x3")]
+
+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
+
+main = print =<< glpSolveVars mipDefaults lp
diff --git a/glpk-hs.cabal b/glpk-hs.cabal
new file mode 100644
--- /dev/null
+++ b/glpk-hs.cabal
@@ -0,0 +1,33 @@
+Name:           glpk-hs
+Version:        0.0.0
+Author:         Louis Wasserman
+License:        GPL
+License-file:   LICENSE
+Maintainer:     Louis Wasserman <wasserman.louis@gmail.com>
+Stability:      experimental
+Synopsis:       Comprehensive GLPK linear programming bindings
+Description:
+    Friendly interface to GLPK's linear programming and mixed integer programming features.  To design a linear programming problem,
+    use "Data.LinearProgram.LPMonad" to construct the constraints and specifications.  Linear functions are essentially specified
+    as @Data.Map@s from variables to their coefficients, and functions for manipulating them are available in "Data.LinFunc".
+    Then "Data.LinearProgram.GLPK" provides facilities for using the GLPK solver system on your problem, with a sizable number
+    of options available.
+
+Category:      Math
+
+cabal-version:  >= 1.2
+build-type:     Simple
+
+extra-source-files: examples/example1.hs
+
+Build-Depends:    base >= 3 && < 5, array, containers, mtl
+Exposed-modules:  Data.LinFunc,
+                  Data.LinearProgram,
+                  Data.LinearProgram.GLPK,
+                  Data.LinearProgram.LPMonad
+Other-modules:    Data.LinearProgram.GLPK.Internal,
+                  Data.LinearProgram.Spec,
+                  Data.LinearProgram.Types,
+                  Data.LinFunc.Class
+c-sources:        glpk/glpk.c
+extra-libraries:  glpk
diff --git a/glpk/glpk.c b/glpk/glpk.c
new file mode 100644
--- /dev/null
+++ b/glpk/glpk.c
@@ -0,0 +1,124 @@
+#include <glpk.h>
+// #include <stdio.h>
+// #include <stdlib.h>
+
+glp_prob *c_glp_create_prob(){
+  	glp_prob *lp;
+	lp = glp_create_prob();
+	return lp;
+}
+
+void c_glp_set_obj_name(glp_prob *lp, const char *name){
+  	glp_set_obj_name(lp, name);
+}
+
+void c_glp_set_obj_dir(glp_prob *lp, int dir){
+  	glp_set_obj_dir(lp, dir ? GLP_MAX : GLP_MIN);
+}
+
+int c_glp_add_rows(glp_prob *lp, int nrows){
+  	return glp_add_rows(lp, nrows);
+}
+
+int c_glp_add_cols(glp_prob *lp, int ncols){
+  	return glp_add_cols(lp, ncols);
+}
+
+void c_glp_set_obj_coef(glp_prob *lp, int j, double coef){
+  	glp_set_obj_coef(lp, j, coef);
+}
+
+void c_glp_set_row_name(glp_prob *lp, int i, const char * name){
+  	glp_set_row_name(lp, i, name);
+}
+
+void c_glp_set_col_name(glp_prob *lp, int i, const char * name){
+	glp_set_col_name(lp, i, name);
+}
+
+void c_glp_set_row_bnds(glp_prob *lp, int i, int type, double lb, double ub){
+	glp_set_row_bnds(lp, i, type, lb, ub);
+}
+
+void c_glp_set_col_bnds(glp_prob *lp, int i, int type, double lb, double ub){
+	glp_set_col_bnds(lp, i, type, lb, ub);
+}
+
+void c_glp_set_mat_row(glp_prob *lp, int i, int len, const int ind[], const double val[]){
+	glp_set_mat_row(lp, i, len, ind, val);
+}
+
+void c_glp_delete_prob(glp_prob *lp){
+  	glp_delete_prob(lp);
+}
+
+void c_glp_create_index(glp_prob *lp){
+  	glp_create_index(lp);
+}
+
+int c_glp_find_row(glp_prob *lp, const char *name){
+  	return glp_find_row(lp, name);
+}
+
+int c_glp_find_col(glp_prob *lp, const char *name){
+  	return glp_find_col(lp, name);
+}
+
+int c_glp_solve_simplex(glp_prob *lp, int msg_lev, int tm_lim, int presolve){
+	glp_smcp smcp;
+	glp_init_smcp (&smcp);
+	smcp.msg_lev = msg_lev;
+	smcp.tm_lim = tm_lim;
+	smcp.presolve = presolve ? GLP_ON : GLP_OFF;
+	glp_adv_basis(lp, 0);
+	return glp_simplex(lp, &smcp);
+}
+
+double c_glp_get_obj_val(glp_prob *lp){
+  	return glp_get_obj_val(lp);
+}
+
+double c_glp_get_row_prim(glp_prob *lp, int i){
+  	return glp_get_row_prim(lp, i);
+}
+
+double c_glp_get_col_prim(glp_prob *lp, int i){
+  	return glp_get_col_prim(lp, i);
+}
+
+void c_glp_set_col_kind(glp_prob *lp, int j, int kind){
+	glp_set_col_kind(lp, j, kind);
+}
+
+int c_glp_mip_solve(glp_prob *lp, int msg_lev, int br_tech, int bt_tech, int pp_tech,
+		     	int fp_heur, int tm_lim, int cuts, double mip_gap, int presolve){
+  	glp_iocp iocp;
+// 	printf ("%d %d %d time\n", msg_lev, br_tech, tm_lim);
+	glp_init_iocp(&iocp);
+	iocp.msg_lev = msg_lev;
+	iocp.br_tech = br_tech;
+	iocp.bt_tech = bt_tech;
+	iocp.pp_tech = pp_tech;
+	iocp.fp_heur = fp_heur;
+	iocp.gmi_cuts = cuts & 1 ? GLP_ON : GLP_OFF;
+	iocp.mir_cuts = cuts & 2 ? GLP_ON : GLP_OFF;
+	iocp.cov_cuts = cuts & 4 ? GLP_ON : GLP_OFF;
+	iocp.clq_cuts = cuts & 8 ? GLP_ON : GLP_OFF;
+	iocp.mip_gap = mip_gap;
+	iocp.tm_lim = tm_lim;
+// 	printf ("%d %d %d time\n", msg_lev, br_tech, tm_lim);
+	iocp.presolve = presolve ? GLP_ON : GLP_OFF;
+	return glp_intopt(lp, &iocp);
+}
+
+double c_glp_mip_obj_val (glp_prob *mip){
+  	return glp_mip_obj_val(mip);
+}
+
+double c_glp_mip_row_val (glp_prob *mip, int i){
+  	return glp_mip_row_val(mip, i);
+}
+
+double c_glp_mip_col_val (glp_prob *mip, int j){
+  	return glp_mip_col_val(mip, j);
+}