diff --git a/Data/LinFunc.hs b/Data/LinFunc.hs
deleted file mode 100644
--- a/Data/LinFunc.hs
+++ /dev/null
@@ -1,121 +0,0 @@
-{-# 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
deleted file mode 100644
--- a/Data/LinFunc/Class.hs
+++ /dev/null
@@ -1,19 +0,0 @@
-{-# 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
--- a/Data/LinearProgram.hs
+++ b/Data/LinearProgram.hs
@@ -1,6 +1,8 @@
-module Data.LinearProgram (module Data.LinearProgram.Spec, module Data.LinearProgram.Types,
-	module Data.LinFunc) where
+module Data.LinearProgram (
+	module Data.LinearProgram.Common,
+	module Data.LinearProgram.LPMonad,
+	module Data.LinearProgram.GLPK) where
 
-import Data.LinearProgram.Spec
-import Data.LinearProgram.Types
-import Data.LinFunc
+import Data.LinearProgram.GLPK
+import Data.LinearProgram.LPMonad
+import Data.LinearProgram.Common
diff --git a/Data/LinearProgram/Common.hs b/Data/LinearProgram/Common.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/Common.hs
@@ -0,0 +1,8 @@
+module Data.LinearProgram.Common (
+	module Data.LinearProgram.Spec,
+	module Data.LinearProgram.LinFunc,
+	module Data.LinearProgram.Types) where
+
+import Data.LinearProgram.Spec
+import Data.LinearProgram.LinFunc
+import Data.LinearProgram.Types
diff --git a/Data/LinearProgram/GLPK.hs b/Data/LinearProgram/GLPK.hs
--- a/Data/LinearProgram/GLPK.hs
+++ b/Data/LinearProgram/GLPK.hs
@@ -1,138 +1,8 @@
-{-# OPTIONS -funbox-strict-fields #-}
-{-# LANGUAGE TupleSections, RecordWildCards #-}
-
-module Data.LinearProgram.GLPK (GLPOpts(..), MsgLev(..), BranchingTechnique(..),
-	BacktrackTechnique(..), Preprocessing(..), Cuts(..), ReturnCode(..), 
-	simplexDefaults, mipDefaults, glpSolveVars, glpSolveAll) where
-
-import Control.Monad
-import Control.Monad.Trans
-
--- import Debug.Trace
-
-import Data.Map
-import Data.Maybe (catMaybes)
-import Data.LinearProgram.Spec
-import Data.LinearProgram.Types
-import Data.LinearProgram.GLPK.Internal
-
-import Data.Time.Clock
--- import System.Time
-
-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, Show 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
-
--- | 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, Show v, Real c) => GLPOpts -> LP v c -> IO (ReturnCode, Maybe (Double, Map v Double, Map String Double))
-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 [maybe (return Nothing) (\ nam -> do
-					val <- getRowPrim i
-					return (Just (nam, val))) nam
-					| (i, Constr nam _ _) <- zip [0..] constraints]
-		return (Just (obj, fromDistinctAscList vals, fromDistinctAscList (catMaybes 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 [maybe (return Nothing) (\ nam -> do
-					val <- mipRowVal i
-					return (Just (nam, val))) nam
-					| (i, Constr nam _ _) <- zip [0..] constraints]
-		return (Just (obj, fromDistinctAscList vals, fromDistinctAscList (catMaybes rows)))) vars
-
-doGLP :: (Ord v, Show 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
-	return (success, guard (gaveAnswer success) >> Just vars)
-doGLP MipOpts{..} lp = do
-	vars <- writeProblem lp
--- 	time <- getTime
--- 	solveSimplex msgLev tmLim presolve
--- 	time' <- getTime
-	let tmLim' = tmLim  --- round (toRational (time' `diffUTCTime` time))
-	success <- mipSolve msgLev brTech btTech ppTech fpHeur cuts mipGap (fromIntegral tmLim') presolve
-	return (success, guard (gaveAnswer success) >> Just vars) --(if success then Just vars else Nothing)
--- 	where	getTime = liftIO getCurrentTime
-
-writeProblem :: (Ord v, Show v, Real c) => LP v c -> GLPK (Map v Int)
-writeProblem LP{..} = do
-	setObjectiveDirection direction
-	i0 <- addCols nVars
-	let allVars' = fmap (i0 +) allVars
-	sequence_ [setColName i (show v) | (v, i) <- assocs allVars']
-	sequence_ [setObjCoef i v | (i, v) <- elems $ intersectionWith (,) allVars' objective]
-	j0 <- addRows (length constraints)
-	sequence_ [do	case lab of
-				Nothing	-> return ()
-				Just n	-> setRowName j n
-			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]
--- 	writeLP
-	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);
-	#-}
+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
new file mode 100644
--- /dev/null
+++ b/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/Data/LinearProgram/GLPK/IO.hs b/Data/LinearProgram/GLPK/IO.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/GLPK/IO.hs
@@ -0,0 +1,19 @@
+module Data.LinearProgram.GLPK.IO where
+
+-- import Control.Monad.Trans
+
+-- import Data.Map
+
+import Data.LinearProgram.Common
+import Data.LinearProgram.LPMonad.Internal
+
+import Data.LinearProgram.GLPK.Common
+import Data.LinearProgram.GLPK.IO.Internal
+
+-- | Read a linear program from a file in CPLEX LP format.
+readLP :: FilePath -> IO (LP String Double)
+readLP = runGLPK . readGLP_LP
+
+-- | 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 . writeGLP_LP file
diff --git a/Data/LinearProgram/GLPK/IO/Internal.hs b/Data/LinearProgram/GLPK/IO/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/GLPK/IO/Internal.hs
@@ -0,0 +1,119 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+module Data.LinearProgram.GLPK.IO.Internal (readGLP_LP, writeGLP_LP) where
+
+import Control.Monad
+import Control.Monad.Trans
+
+
+import Data.Map hiding (map)
+
+import Data.LinearProgram.Common
+import Data.LinearProgram.GLPK.Common
+import Data.LinearProgram.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_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)
+
+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 (lb i)
+		4	-> liftM2 Bound (lb i) (ub i)
+		5	-> 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
+		coefs <- liftIO $ mallocArray m
+		sequence [do
+			k <- liftM fromIntegral $ GLP $ \ lp -> glpGetMatRow lp (fromIntegral i) ixs coefs
+			ixsL <- liftIO $ peekArray k ixs
+			coefsL <- liftIO $ peekArray k coefs
+			return (i, zip (map fromIntegral ixsL) (map realToFrac coefsL))
+			| i <- [1..n]]
+
+readGLP_LP :: FilePath -> GLPK (LP String Double)
+readGLP_LP 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]
+
+writeGLP_LP :: (Show v, Ord v, Real c) => FilePath -> LP v c -> GLPK ()
+writeGLP_LP file lp = do
+	writeProblem lp
+	writeLP file
diff --git a/Data/LinearProgram/GLPK/Internal.hs b/Data/LinearProgram/GLPK/Internal.hs
--- a/Data/LinearProgram/GLPK/Internal.hs
+++ b/Data/LinearProgram/GLPK/Internal.hs
@@ -1,36 +1,21 @@
-{-# LANGUAGE ScopedTypeVariables, EmptyDataDecls, ForeignFunctionInterface #-}
-module Data.LinearProgram.GLPK.Internal (GLPK, MsgLev (..), Preprocessing (..), Direction(..), BacktrackTechnique(..),
-	BranchingTechnique(..), Cuts(..), ReturnCode(..), gaveAnswer, runGLPK, writeLP, addCols,
+{-# LANGUAGE RecordWildCards, ScopedTypeVariables, EmptyDataDecls, ForeignFunctionInterface #-}
+module Data.LinearProgram.GLPK.Internal (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 Control.Monad.Trans
 
--- import Debug.Trace
-
 import Foreign.Ptr
 import Foreign.C
-import Foreign.ForeignPtr
 import Foreign.Marshal.Array
 
 import Data.Bits
+import Data.Map hiding (map)
 -- import Data.Bounds
-import Data.LinearProgram.Types
-
-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]
+import Data.LinearProgram.Common
+import Data.LinearProgram.GLPK.Types
 
-foreign import ccall unsafe "c_glp_create_prob" glpCreateProb :: IO (Ptr GlpProb)
 -- 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_add_rows" glpAddRows :: Ptr GlpProb -> CInt -> IO CInt
@@ -41,7 +26,6 @@
 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_delete_prob" glpDelProb :: FunPtr (Ptr GlpProb -> 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
@@ -55,31 +39,10 @@
 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_write_lp" glpWriteLP :: Ptr GlpProb -> IO ()
 
-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)
-
-writeLP :: GLPK ()
-writeLP = GLP glpWriteLP
-
-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)
+	(fromIntegral $ 1 + fromEnum dir)
 
 addRows :: Int -> GLPK Int
 addRows n = GLP $ liftM fromIntegral . flip glpAddRows (fromIntegral n)
@@ -127,8 +90,6 @@
 findCol :: String -> GLPK Int
 findCol nam = GLP $ liftM fromIntegral . withCString nam . glpFindCol
 
-data MsgLev = MsgOff | MsgErr | MsgOn | MsgAll
-
 solveSimplex :: MsgLev -> Int -> Bool -> GLPK ReturnCode
 solveSimplex msglev tmLim presolve = GLP $ \ lp -> liftM (toEnum . fromIntegral) $ glpSolveSimplex lp
 	(getMsgLev msglev)
@@ -137,11 +98,7 @@
 	where	tmLim' = fromIntegral (tmLim * 1000)
 
 getMsgLev :: MsgLev -> CInt
-getMsgLev msglev = case msglev of
-	MsgOff	-> 0
-	MsgErr	-> 1
-	MsgOn	-> 2
-	MsgAll	-> 3
+getMsgLev = fromIntegral . fromEnum
 
 getObjVal :: GLPK Double
 getObjVal = liftM realToFrac $ GLP glpGetObjVal
@@ -153,15 +110,7 @@
 getColPrim i = liftM realToFrac $ GLP (`glpGetColPrim` fromIntegral i)
 
 setColKind :: Int -> VarKind -> GLPK ()
-setColKind i kind = GLP $ \ lp -> glpSetColKind lp (fromIntegral i) (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)
+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
@@ -201,3 +150,30 @@
 
 mipColVal :: Int -> GLPK Double
 mipColVal i = liftM realToFrac $ GLP (`glpMIPColVal` fromIntegral i)
+
+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	case lab of
+				Nothing	-> return ()
+				Just n	-> setRowName j n
+			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]
+-- 	writeLP
+	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
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/GLPK/Solver.hs
@@ -0,0 +1,121 @@
+{-# OPTIONS -funbox-strict-fields #-}
+{-# LANGUAGE TupleSections, RecordWildCards #-}
+
+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 Control.Monad.Trans
+
+-- import Debug.Trace
+
+import Data.Map
+import Data.Maybe (catMaybes)
+import Data.LinearProgram.Common
+import Data.LinearProgram.GLPK.Internal
+import Data.LinearProgram.GLPK.Types
+
+import Data.Time.Clock
+-- import System.Time
+
+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}
+
+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
+
+-- | 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
+
+-- | 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) (getRowPrim i)
+					| (i, c) <- zip [1..] constraints]
+		return (Just (obj, fromDistinctAscList vals, rows))) vars
+
+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
+	return (success, guard (gaveAnswer success) >> Just vars)
+doGLP MipOpts{..} lp = do
+	vars <- writeProblem lp
+-- 	time <- getTime
+-- 	solveSimplex msgLev tmLim presolve
+-- 	time' <- getTime
+	let tmLim' = tmLim  --- round (toRational (time' `diffUTCTime` time))
+	success <- mipSolve msgLev brTech btTech ppTech fpHeur cuts mipGap (fromIntegral tmLim') presolve
+	return (success, guard (gaveAnswer success) >> Just vars) --(if success then Just vars else Nothing)
+-- 	where	getTime = liftIO getCurrentTime
+		
+{-# 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/Types.hs b/Data/LinearProgram/GLPK/Types.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/GLPK/Types.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE EmptyDataDecls, ForeignFunctionInterface #-}
+
+module Data.LinearProgram.GLPK.Types where
+
+import Control.Monad.Trans
+
+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 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/LPMonad.hs b/Data/LinearProgram/LPMonad.hs
--- a/Data/LinearProgram/LPMonad.hs
+++ b/Data/LinearProgram/LPMonad.hs
@@ -1,88 +1,75 @@
-{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE FlexibleContexts, RecordWildCards #-}
 
-module Data.LinearProgram.LPMonad where
+-- | 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.
+module Data.LinearProgram.LPMonad (
+	module Data.LinearProgram.LPMonad.Internal,
+	-- * Solvers
+	quickSolveMIP,
+	quickSolveLP,
+	glpSolve,
+	quickSolveMIP',
+	quickSolveLP',
+	glpSolve') where
 
 import Control.Monad.State.Strict
+import Control.Monad.Identity
 
 import Data.Map
 import Data.Monoid
--- import Data.Bounds
 
-import Data.LinFunc
-import Data.LinearProgram.Types
-import Data.LinearProgram.Spec
+import Data.LinearProgram.Common
+import Data.LinearProgram.LPMonad.Internal
 
--- | A 'State' monad used for the construction of a linear program.
-type LPM v c = State (LP v c)
+import Data.LinearProgram.GLPK.Solver
+import Data.LinearProgram.GLPK.IO
 
--- | 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}
+{-# 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
 
-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}
+{-# 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
 
-setObjective :: LinFunc v c -> LPM v c ()
-setObjective obj = modify setObj where
-	setObj lp = lp{objective = obj}
+{-# 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
 
-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)
+{-# 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
 
-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}
+{-# SPECIALIZE writeLPToFile :: (Ord v, Show v, Real c) => FilePath -> LPT v c IO () #-}
+writeLPToFile :: (Ord v, Show v, Real c, MonadState (LP v c) m, MonadIO m) =>
+	FilePath -> m ()
+writeLPToFile file = get >>= liftIO . writeLP file 
 
-setVarKind :: Ord v => v -> VarKind -> LPM v c ()
-setVarKind v k = modify setK where
-	setK lp@LP{..} = lp{varTypes = insertWith mappend v k varTypes}
+{-# SPECIALIZE readLPFromFile :: FilePath -> LPT String Double IO () #-}
+readLPFromFile :: (MonadState (LP String Double) m, MonadIO m) =>
+	FilePath -> m ()
+readLPFromFile file = put =<< liftIO (readLP file)
diff --git a/Data/LinearProgram/LPMonad/Internal.hs b/Data/LinearProgram/LPMonad/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/LPMonad/Internal.hs
@@ -0,0 +1,200 @@
+{-# LANGUAGE FlexibleContexts, RecordWildCards #-}
+
+module Data.LinearProgram.LPMonad.Internal (
+-- 	module Data.LinearProgram.Common,
+	-- * Monad definitions
+	LPM,
+	LPT,
+	runLPM,
+	runLPT,
+	execLPM,
+	execLPT,
+	evalLPM,
+	evalLPT,
+	-- * 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) where
+
+import Control.Monad.State.Strict
+import Control.Monad.Identity
+
+import Data.Map
+import Data.Monoid
+
+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, Module r c) => LPM v c a -> (a, LP v c)
+runLPM = runIdentity . runLPT
+
+runLPT :: (Ord v, Module r 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, Module r c) => LPM v c a -> LP v c
+execLPM = runIdentity . execLPT
+
+-- | Constructs a linear programming problem in the specified monad.
+execLPT :: (Ord v, Module r 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, Module r c) => LPM v c a -> a
+evalLPM = runIdentity . evalLPT
+
+-- | Runs the specified operation in the linear programming monad transformer.
+evalLPT :: (Monad m, Ord v, Module r c) => 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, Module r c) => LinFunc v c -> LinFunc v c -> LPM v c (),
+	(Ord v, Module r c, Monad m) => LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+{-# SPECIALIZE leq :: (Ord v, Module r c) => LinFunc v c -> LinFunc v c -> LPM v c (),
+	(Ord v, Module r c, Monad m) => LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+{-# SPECIALIZE geq :: (Ord v, Module r c) => LinFunc v c -> LinFunc v c -> LPM v c (),
+	(Ord v, Module r c, Monad m) => LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+-- | Specifies the relationship between two functions in the variables.
+equal, leq, geq :: (Ord v, Module r 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, Module r c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
+	(Ord v, Module r c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+{-# SPECIALIZE geq' :: (Ord v, Module r c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
+	(Ord v, Module r c, Monad m) => String -> LinFunc v c -> LinFunc v c -> LPT v c m () #-}
+{-# SPECIALIZE leq' :: (Ord v, Module r c) => String -> LinFunc v c -> LinFunc v c -> LPM v c (),
+	(Ord v, Module r 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, Module r 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 labelled 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 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.
+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.
+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, Module r c) => LinFunc v c -> LPM v c (),
+	(Ord v, Module r c, Monad m) => LinFunc v c -> LPT v c m () #-}
+-- | Adds this function to the objective function.
+addObjective :: (Ord v, Module r 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.
+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/Data/LinearProgram/LinFunc.hs b/Data/LinearProgram/LinFunc.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/LinFunc.hs
@@ -0,0 +1,122 @@
+{-# LANGUAGE UndecidableInstances, FlexibleInstances, MultiParamTypeClasses #-}
+
+module Data.LinearProgram.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.LinearProgram.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]
+
+{-# 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]
diff --git a/Data/LinearProgram/LinFunc/Class.hs b/Data/LinearProgram/LinFunc/Class.hs
new file mode 100644
--- /dev/null
+++ b/Data/LinearProgram/LinFunc/Class.hs
@@ -0,0 +1,19 @@
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}
+
+module Data.LinearProgram.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/Spec.hs b/Data/LinearProgram/Spec.hs
--- a/Data/LinearProgram/Spec.hs
+++ b/Data/LinearProgram/Spec.hs
@@ -3,7 +3,7 @@
 
 -- import Control.DeepSeq
 -- import Data.Bounds
-import Data.LinFunc
+import Data.LinearProgram.LinFunc
 import Data.LinearProgram.Types
 import Data.Map
 
diff --git a/Data/LinearProgram/Types.hs b/Data/LinearProgram/Types.hs
--- a/Data/LinearProgram/Types.hs
+++ b/Data/LinearProgram/Types.hs
@@ -4,7 +4,7 @@
 
 import Data.Monoid
 
-data VarKind = ContVar | IntVar | BinVar deriving (Eq, Ord, Show, Read)
+data VarKind = ContVar | IntVar | BinVar deriving (Eq, Ord, Enum, Show, Read)
 
 -- instance NFData VarKind
 
@@ -12,7 +12,7 @@
 	mempty = ContVar
 	mappend = max
 
-data Direction = Min | Max deriving (Eq, Ord, Show, Read)
+data Direction = Min | Max deriving (Eq, Ord, Enum, Show, Read)
 
 -- instance NFData Direction
 
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.0.2
+Version:        0.0.3
 Author:         Louis Wasserman
 License:        GPL
 License-file:   LICENSE
@@ -21,13 +21,20 @@
 extra-source-files: examples/example1.hs
 
 Build-Depends:    base >= 3 && < 5, array, containers, mtl, time
-Exposed-modules:  Data.LinFunc,
-                  Data.LinearProgram,
+Exposed-modules:  Data.LinearProgram,
+                  Data.LinearProgram.Common,
                   Data.LinearProgram.GLPK,
+                  Data.LinearProgram.GLPK.Solver,
+                  Data.LinearProgram.GLPK.IO,
                   Data.LinearProgram.LPMonad
-Other-modules:    Data.LinearProgram.GLPK.Internal,
+Other-modules:    Data.LinearProgram.GLPK.Internal,	
+                  Data.LinearProgram.GLPK.Types,
+                  Data.LinearProgram.GLPK.Common,
+                  Data.LinearProgram.GLPK.IO.Internal,
+                  Data.LinearProgram.LPMonad.Internal,
                   Data.LinearProgram.Spec,
                   Data.LinearProgram.Types,
-                  Data.LinFunc.Class
+                  Data.LinearProgram.LinFunc,
+                  Data.LinearProgram.LinFunc.Class
 c-sources:        glpk/glpk.c
 extra-libraries:  glpk
diff --git a/glpk/glpk.c b/glpk/glpk.c
--- a/glpk/glpk.c
+++ b/glpk/glpk.c
@@ -8,10 +8,6 @@
 	return lp;
 }
 
-void c_glp_write_lp (glp_prob *lp){
-  	glp_write_lp(lp, 0, "test.cplex");
-}
-
 void c_glp_set_obj_name(glp_prob *lp, const char *name){
   	glp_set_obj_name(lp, name);
 }
@@ -126,4 +122,72 @@
 
 double c_glp_mip_col_val (glp_prob *mip, int j){
   	return glp_mip_col_val(mip, j);
+}
+
+int c_glp_get_obj_dir (glp_prob *lp){
+  	return glp_get_obj_dir (lp);
+}
+
+int c_glp_get_num_rows (glp_prob *lp){
+  	return glp_get_num_rows (lp);
+}
+
+int c_glp_get_num_cols (glp_prob *lp){
+ 	return glp_get_num_cols (lp);
+}
+
+const char *c_glp_get_row_name (glp_prob *lp, int i){
+  	return glp_get_row_name (lp, i);
+}
+
+const char *c_glp_get_col_name (glp_prob *lp, int i){
+  	return glp_get_col_name (lp, i);
+}
+
+int c_glp_get_col_kind (glp_prob *lp, int i){
+  	return glp_get_col_kind (lp, i);
+}
+
+// double c_dbl_max(){
+//   	return DBL_MAX;
+// }
+
+double c_glp_get_row_lb (glp_prob *lp, int i){
+	return glp_get_row_lb (lp, i);
+}
+
+double c_glp_get_row_ub (glp_prob *lp, int i){
+	return glp_get_row_ub (lp, i);
+}
+
+double c_glp_get_col_lb (glp_prob *lp, int i){
+	return glp_get_col_lb (lp, i);
+}
+
+double c_glp_get_col_ub (glp_prob *lp, int i){
+	return glp_get_col_ub (lp, i);
+}
+
+int c_glp_get_row_type (glp_prob *lp, int i){
+  	return glp_get_row_type(lp, i);
+}
+
+int c_glp_get_col_type (glp_prob *lp, int i){
+  	return glp_get_col_type(lp, i);
+}
+
+double c_glp_get_obj_coef (glp_prob *lp, int j){
+  	return glp_get_obj_coef(lp, j);
+}
+
+int c_glp_get_mat_row (glp_prob *lp, int i, int ind[], double val[]){
+  	return c_glp_get_mat_row (lp, i, ind, val);
+}
+
+int c_glp_read_lp(glp_prob *lp, const char * fname){
+  	return glp_read_lp(lp, NULL, fname);
+}
+
+int c_glp_write_lp(glp_prob *lp, const char * fname){
+  	return glp_write_lp(lp, NULL, fname);
 }