diff --git a/Data/Algebra.hs b/Data/Algebra.hs
new file mode 100644
--- /dev/null
+++ b/Data/Algebra.hs
@@ -0,0 +1,32 @@
+-- | 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(..),
+	Module(..),
+	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.Module
diff --git a/Data/Algebra/Group.hs b/Data/Algebra/Group.hs
new file mode 100644
--- /dev/null
+++ b/Data/Algebra/Group.hs
@@ -0,0 +1,102 @@
+{-# 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
new file mode 100644
--- /dev/null
+++ b/Data/Algebra/Module.hs
@@ -0,0 +1,109 @@
+{-# 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, 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, 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
new file mode 100644
--- /dev/null
+++ b/Data/Algebra/Ring.hs
@@ -0,0 +1,58 @@
+{-# 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
+	(*#) = (*)
+
+instance Ring r => Ring (Poly r) where
+	one = [one]
+	(p:ps) *# (q:qs) = (p *# q):(ps *# (q:qs) ^+^ map (p *#) qs)
+	_ *# _ = []
+
+instance Ring r => Ring (a -> r) where
+	one = const one
+	(*#) = liftA2 (*#)
+
+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/Common.hs b/Data/LinearProgram/Common.hs
--- a/Data/LinearProgram/Common.hs
+++ b/Data/LinearProgram/Common.hs
@@ -1,10 +1,12 @@
+-- | 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.LinearProgram.LinFunc,
+	module Data.Algebra,
 	module Data.LinearProgram.Types) where
 
 import Data.LinearProgram.Spec
-import Data.LinearProgram.LinFunc
+import Data.Algebra
 import Data.LinearProgram.Types
 
 import Data.Map
diff --git a/Data/LinearProgram/GLPK/IO.hs b/Data/LinearProgram/GLPK/IO.hs
--- a/Data/LinearProgram/GLPK/IO.hs
+++ b/Data/LinearProgram/GLPK/IO.hs
@@ -1,18 +1,20 @@
+-- | Bindings to the file I/O functions from GLPK, on the CPLEX LP file format.
 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
 
+{-# 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, Ord c, Fractional c, Group c) => FilePath -> IO (LP v c)
+readLP = fmap (mapVars read . mapVals realToFrac) . readLP'
+
 -- | Read a linear program from a file in CPLEX LP format.
-readLP :: FilePath -> IO (LP String Double)
-readLP = runGLPK . readGLPLP
+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 ()
diff --git a/Data/LinearProgram/GLPK/IO/Internal.hs b/Data/LinearProgram/GLPK/IO/Internal.hs
--- a/Data/LinearProgram/GLPK/IO/Internal.hs
+++ b/Data/LinearProgram/GLPK/IO/Internal.hs
@@ -123,7 +123,6 @@
 		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
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,6 +1,6 @@
 {-# LANGUAGE RecordWildCards, ScopedTypeVariables, ForeignFunctionInterface #-}
 module Data.LinearProgram.GLPK.Internal (writeProblem, solveSimplex, mipSolve,
-	getObjVal, getRowPrim, getColPrim, mipObjVal, mipRowVal, mipColVal) where
+	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,
@@ -19,7 +19,9 @@
 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_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 ()
@@ -39,11 +41,17 @@
 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 $ flip glpSetObjDir 
-	(fromIntegral $ 1 + fromEnum dir)
+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)
 
@@ -64,9 +72,11 @@
 	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
@@ -133,6 +143,11 @@
 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
@@ -140,17 +155,17 @@
 	let allVars' = fmap (i0 +) allVars
 	sequence_ [setObjCoef i v | (i, v) <- elems $ intersectionWith (,) allVars' objective]
 	j0 <- addRows (length constraints)
-	sequence_ [do	setMatRow j
+	sequence_ [do	maybe (return ()) (setRowName j) lab
+			setMatRow j
 				[(i, v) | (i, v) <- elems (intersectionWith (,) allVars' f)]
 			setRowBounds j bnds
-				| (j, Constr _ f bnds) <- zip [j0..] constraints]
+				| (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`
diff --git a/Data/LinearProgram/GLPK/Solver.hs b/Data/LinearProgram/GLPK/Solver.hs
--- a/Data/LinearProgram/GLPK/Solver.hs
+++ b/Data/LinearProgram/GLPK/Solver.hs
@@ -1,6 +1,16 @@
 {-# 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(..),
@@ -40,6 +50,8 @@
 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))
@@ -62,6 +74,9 @@
 				| (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.
@@ -89,17 +104,16 @@
 					| (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
-	return (success, guard (gaveAnswer success) >> Just vars)
+	bad <- getBadRay
+	maybe (return (success, guard (gaveAnswer success) >> Just vars)) (fail . show) bad
 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
+	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/LPMonad.hs b/Data/LinearProgram/LPMonad.hs
--- a/Data/LinearProgram/LPMonad.hs
+++ b/Data/LinearProgram/LPMonad.hs
@@ -77,14 +77,16 @@
 {-# SPECIALIZE readLPFromFile :: (Ord v, Read v, Ord c, Fractional c, Module r 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, Ord c, Fractional c, Module r c, MonadState (LP v c) m, MonadIO m) =>
+readLPFromFile :: (Ord v, Read v, Ord c, Fractional c, Group c, MonadState (LP v c) m, MonadIO m) =>
 	FilePath -> m ()
-readLPFromFile = put . mapVars read . mapVals realToFrac <=< liftIO . execLPT . readLPFromFile'
+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' file = put =<< liftIO (readLP file)
+readLPFromFile' = put <=< liftIO . readLP'
diff --git a/Data/LinearProgram/LPMonad/Internal.hs b/Data/LinearProgram/LPMonad/Internal.hs
--- a/Data/LinearProgram/LPMonad/Internal.hs
+++ b/Data/LinearProgram/LPMonad/Internal.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleContexts, RecordWildCards #-}
+{-# LANGUAGE BangPatterns, FlexibleContexts, RecordWildCards #-}
 
 module Data.LinearProgram.LPMonad.Internal (
 -- 	module Data.LinearProgram.Common,
@@ -39,7 +39,8 @@
 	varGeq,
 	varBds,
 	setVarBounds,
-	setVarKind) where
+	setVarKind,
+	newVariables) where
 
 import Control.Monad.State.Strict
 import Control.Monad.Identity
@@ -56,26 +57,26 @@
 -- | 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 :: (Ord v, Group 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 :: (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, Module r c) => LPM v c a -> LP v c
+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, Module r c, Monad m) => LPT v c m a -> m (LP v c)
+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, Module r c) => LPM v c a -> a
+evalLPM :: (Ord v, Group 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 :: (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.
@@ -83,26 +84,26 @@
 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 () #-}
+{-# 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.
-equal, leq, geq :: (Ord v, Module r c, MonadState (LP v c) m) => LinFunc v c -> LinFunc v c -> m ()
+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, 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 () #-}
+{-# 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, Module r c, MonadState (LP v c) m) => String -> LinFunc v c -> LinFunc v c -> m ()
+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'
@@ -128,6 +129,22 @@
 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.
+newVariables :: (MonadState (LP v c) m, Ord v, Enum v) => Int -> m [v]
+newVariables !k = do	LP{..} <- get
+			let allVars0 = fmap (const ()) objective `union`
+				unions [fmap (const ()) f | Constr _ f _ <- constraints] `union`
+				fmap (const ()) varBounds `union` fmap (const ()) varTypes
+			case minViewWithKey allVars0 of
+				Nothing	-> return $ take k [toEnum 0..]
+				Just ((start, _), _)
+					-> return $ take k $ tail [start..]
+
 {-# SPECIALIZE 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 (),
@@ -173,10 +190,10 @@
 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 () #-}
+{-# 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, Module r c, MonadState (LP v c) m) => LinFunc v c -> m ()
+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}
 		
diff --git a/Data/LinearProgram/LinFunc.hs b/Data/LinearProgram/LinFunc.hs
deleted file mode 100644
--- a/Data/LinearProgram/LinFunc.hs
+++ /dev/null
@@ -1,120 +0,0 @@
-{-# LANGUAGE UndecidableInstances, FlexibleInstances, MultiParamTypeClasses #-}
-
-module Data.LinearProgram.LinFunc (LinFunc, Module(..), var, varSum, (*&), vsum, combination, linCombination) where
-
-import Control.Applicative
-
-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.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
-	(^+^) = liftA2 (^+^)
-	(^-^) = liftA2 (^-^)
-	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
deleted file mode 100644
--- a/Data/LinearProgram/LinFunc/Class.hs
+++ /dev/null
@@ -1,19 +0,0 @@
-{-# 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
@@ -11,16 +11,23 @@
 
 import Text.ParserCombinators.ReadP
 
-import Data.LinearProgram.LinFunc
+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)
 
@@ -110,14 +117,17 @@
 -- | 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, Module r c) => (v -> v') -> LP v c -> LP v' c
+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],
@@ -126,7 +136,7 @@
 
 -- | Applies the specified function to the constants in the linear program.  This is only safe
 -- for a monotonic function.
-mapVals :: (Ord c', Module r c') => (c -> c') -> LP v c -> LP v c'
+mapVals :: (c -> c') -> LP v c -> LP v c'
 mapVals = fmap
 
 -- instance (NFData v, NFData c) => NFData (Constraint v c) where
diff --git a/Data/LinearProgram/Types.hs b/Data/LinearProgram/Types.hs
--- a/Data/LinearProgram/Types.hs
+++ b/Data/LinearProgram/Types.hs
@@ -29,7 +29,22 @@
 	bd `mappend` Free = bd
 	Equ a `mappend` Equ b
 		| a == b	= Equ a
-		| otherwise	= infeasible
+	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
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.1.0
+Version:        0.2.0
 Author:         Louis Wasserman
 License:        GPL
 License-file:   LICENSE
@@ -7,7 +7,10 @@
 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,
+    Friendly interface to GLPK's linear programming and mixed integer programming features.  Intended for easy extensibility,
+    with a general, pure-Haskell representation of linear programs.  Also includes usefully general algebraic structures.
+    
+    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
@@ -20,13 +23,14 @@
 
 extra-source-files: examples/example1.hs
 
-Build-Depends:    base >= 4 && < 5, array, containers, mtl, time
+Build-Depends:    base >= 4 && < 5, array, containers, mtl
 Exposed-modules:  Data.LinearProgram,
                   Data.LinearProgram.Common,
                   Data.LinearProgram.GLPK,
                   Data.LinearProgram.GLPK.Solver,
                   Data.LinearProgram.GLPK.IO,
-                  Data.LinearProgram.LPMonad
+                  Data.LinearProgram.LPMonad,
+                  Data.Algebra
 Other-modules:    Data.LinearProgram.GLPK.Internal,	
                   Data.LinearProgram.GLPK.Types,
                   Data.LinearProgram.GLPK.Common,
@@ -34,7 +38,8 @@
                   Data.LinearProgram.LPMonad.Internal,
                   Data.LinearProgram.Spec,
                   Data.LinearProgram.Types,
-                  Data.LinearProgram.LinFunc,
-                  Data.LinearProgram.LinFunc.Class
+                  Data.Algebra.Group,
+                  Data.Algebra.Ring,
+                  Data.Algebra.Module
 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
@@ -12,8 +12,20 @@
   	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);
+void c_glp_maximize(glp_prob *lp){
+  	glp_set_obj_dir(lp, GLP_MAX);
+}
+
+void c_glp_minimize(glp_prob *lp){
+  	glp_set_obj_dir(lp, GLP_MIN);
+}
+
+// void c_glp_set_obj_dir(glp_prob *lp, int dir){
+//   	glp_set_obj_dir(lp, dir);
+// }
+
+int c_glp_get_bad_ray(glp_prob *lp){
+  	return glp_get_unbnd_ray(lp);
 }
 
 int c_glp_add_rows(glp_prob *lp, int nrows){