diff --git a/Makefile b/Makefile
new file mode 100644
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,11 @@
+run-test:	update-test
+	runhaskell Setup configure --user --enable-tests
+	runhaskell Setup build
+	runhaskell Setup haddock
+	runhaskell Setup test comfort-glpk-test --show-details=streaming
+
+	runhaskell Setup configure --user -fdebug
+	runhaskell Setup build
+
+update-test:
+	doctest-extract-0.1 -i src/ -o test/ --executable-main=Main.hs Numeric.GLPK Numeric.GLPK.Monad
diff --git a/comfort-glpk.cabal b/comfort-glpk.cabal
--- a/comfort-glpk.cabal
+++ b/comfort-glpk.cabal
@@ -1,6 +1,6 @@
 Cabal-Version:    2.2
 Name:             comfort-glpk
-Version:          0.0.1
+Version:          0.1
 License:          BSD-3-Clause
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -32,8 +32,16 @@
   .
   Alternatives: @coinor-clp@, @hmatrix-glpk@, @glpk-hs@
 
+Extra-Source-Files:
+  Makefile
+
+Flag debug
+  Description: Enable debug output
+  Default:     False
+  Manual:      True
+
 Source-Repository this
-  Tag:         0.0.1
+  Tag:         0.1
   Type:        darcs
   Location:    https://hub.darcs.net/thielema/comfort-glpk/
 
@@ -44,6 +52,7 @@
 Library
   Build-Depends:
     glpk-headers >=0.4.1 && <0.6,
+    linear-programming >=0.0 && <0.1,
     comfort-array >=0.4 && <0.6,
     deepseq >=1.3 && <1.5,
     transformers >=0.3 && <0.7,
@@ -52,20 +61,27 @@
     base >=4.5 && <5
 
   GHC-Options:      -Wall
-  Hs-Source-Dirs:   src
   Default-Language: Haskell98
+  Hs-Source-Dirs:   src
+  If flag(debug)
+    Hs-Source-Dirs: src/debug-on
+  Else
+    Hs-Source-Dirs: src/debug-off
   Extra-Libraries: glpk
   Exposed-modules:
     Numeric.GLPK
     Numeric.GLPK.Monad
   Other-Modules:
     Numeric.GLPK.Private
+    Numeric.GLPK.Debug
 
 Test-Suite comfort-glpk-test
   Type:             exitcode-stdio-1.0
   Build-Depends:
     comfort-glpk,
+    linear-programming,
     comfort-array >=0.5.2,
+    transformers,
     non-empty,
     utility-ht >=0.0.17,
     doctest-exitcode-stdio >=0.0 && <0.1,
@@ -79,6 +95,6 @@
   Default-Language: Haskell98
   Main-Is: Main.hs
   Other-Modules:
-    Test.Numeric.GLPK.Generator
+    Test.Numeric.GLPK.Utility
     Test.Numeric.GLPK.Monad
     Test.Numeric.GLPK
diff --git a/src/Numeric/GLPK.hs b/src/Numeric/GLPK.hs
--- a/src/Numeric/GLPK.hs
+++ b/src/Numeric/GLPK.hs
@@ -23,53 +23,38 @@
 instead of @[i >=. a, i <=. b]@ use @i >=<. (a,b)@.
 -}
 module Numeric.GLPK (
-   Term(..),
+   Term,
    Bound(..),
    Inequality(..),
-   free, (<=.), (>=.), (==.), (>=<.),
-   NoSolutionType(..),
+   LP.free, (LP.<=.), (LP.>=.), (LP.==.), (LP.>=<.),
+   FailureType(..),
    SolutionType(..),
-   Solution,
+   Result,
    Constraints,
    Direction(..),
    Objective,
    Bounds,
    (.*),
-   objectiveFromTerms,
+   LP.objectiveFromTerms,
    simplex,
-   simplexMulti,
-   simplexSuccessive,
    exact,
-   exactMulti,
-   exactSuccessive,
    interior,
-   interiorMulti,
-   interiorSuccessive,
-
-   solveSuccessive,
-
-   FormatIdentifier,
-   formatMathProg,
    ) where
 
 import qualified Math.Programming.Glpk.Header as FFI
+import qualified Numeric.GLPK.Debug as Debug
+import qualified Numeric.LinearProgramming.Common as LP
 import Numeric.GLPK.Private
+import Numeric.LinearProgramming.Common
+         (Bound(..), Inequality(Inequality),
+          Bounds, Direction(..), Objective, (.*))
 
 import qualified Data.Array.Comfort.Storable.Mutable as Mutable
 import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import qualified Data.NonEmpty as NonEmpty
-import qualified Data.List as List
-import Data.Array.Comfort.Storable (Array)
-import Data.Tuple.HT (mapFst, mapSnd)
-import Data.Traversable (for)
 import Data.Foldable (for_)
 
-import Text.Printf (printf)
-
-import qualified Control.Monad.Trans.Except as ME
-import qualified Control.Monad.Trans.State as MS
-import Control.Monad (void, when)
+import Control.Monad (void)
 import Control.Applicative (liftA2)
 import Control.Exception (bracket)
 
@@ -80,7 +65,7 @@
 
 
 {- $setup
->>> import qualified Test.Numeric.GLPK.Generator as TestLP
+>>> import qualified Numeric.LinearProgramming.Test as TestLP
 >>> import qualified Numeric.GLPK as LP
 >>> import Numeric.GLPK ((.*), (<=.), (==.))
 >>>
@@ -107,176 +92,37 @@
 -}
 
 
-infix 7 .*
-
-(.*) :: Double -> ix -> Term ix
-(.*) = Term
-
-
-infix 4 <=., >=., >=<., ==.
-
-(<=.), (>=.), (==.) :: x -> Double -> Inequality x
-x <=. bnd = Inequality x $ LessEqual bnd
-x >=. bnd = Inequality x $ GreaterEqual bnd
-x ==. bnd = Inequality x $ Equal bnd
-
-(>=<.) :: x -> (Double,Double) -> Inequality x
-x >=<. bnd = Inequality x $ uncurry Between bnd
-
-free :: x -> Inequality x
-free x = Inequality x Free
-
-
-
-objectiveFromTerms ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix) => sh -> [Term ix] -> Objective sh
-objectiveFromTerms sh =
-   Array.fromAssociations 0 sh . map (\(Term x ix) -> (ix,x))
-
-
 {- |
 >>> case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
 Right (Optimal,(28.0,(5.0,0.0,4.0)))
 
 prop> \target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (1,1) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (TestLP.approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False
 prop> \(QC.Positive posWeight) (QC.Positive negWeight) target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (posWeight,negWeight) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, _) -> True; _ -> False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, _) -> True; _ -> False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.1 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, (opt,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> let val = TestLP.scalarProduct obj $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin) in (case dir of LP.Minimize -> opt-0.01 <= val; LP.Maximize -> opt+0.01 >= val); _ -> QC.property False
 -}
 simplex ::
    (Shape.Indexed sh, Shape.Index sh ~ ix) =>
    Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) -> Solution sh
+   (Direction, Objective sh) -> Result sh
 simplex = solve (flip FFI.glp_simplex nullPtr)
 
-{-# DEPRECATED simplexMulti "use GLPK.Monad instead" #-}
-{-# DEPRECATED exactMulti "use GLPK.Monad instead" #-}
-{-# DEPRECATED interiorMulti "run 'interior' in Either monad instead" #-}
-{- |
-Optimize for one objective after another.
-That is, if the first optimization succeeds
-then the optimum is fixed as constraint
-and the optimization is continued with respect to the second objective and so on.
-The iteration fails if one optimization fails.
-The obtained objective values are returned as well.
-Their number equals the number of attempted optimizations.
 
-The last objective value is included in the Solution value.
-This is a bit inconsistent,
-but this way you have a warranty that there is an objective value
-if the optimization is successful.
-
-The objectives are expected as 'Term's
-because after successful optimization step
-they are used as (sparse) constraints.
-It's also easy to assert that the same array shape
-is used for all objectives.
-
-The function does not work reliably,
-because an added objective can make the system infeasible
-due to rounding errors.
-E.g. a non-negative objective can become very small but negative.
-
-
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case LP.simplexMulti bounds constrs (Array.shape origin) objs of (_, Right (LP.Optimal, _)) -> QC.property True; result -> QC.counterexample (show result) False
-
-The same property fails for 'exactMulti' and 'interiorMulti'.
-I guess, due to rounding errors.
--}
-simplexMulti, exactMulti, interiorMulti ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
-   Bounds ix -> Constraints ix ->
-   sh -> NonEmpty.T [] (Direction, [Term ix]) -> ([Double], Solution sh)
-simplexMulti = solveMulti . simplex
-exactMulti = solveMulti . exact
-interiorMulti = solveMulti . interior
-
-solveMulti ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
-   (Constraints ix -> (Direction, Objective sh) -> Solution sh) ->
-   Constraints ix ->
-   sh -> NonEmpty.T [] (Direction, [Term ix]) -> ([Double], Solution sh)
-solveMulti solver constrs0 sh (NonEmpty.Cons obj0 objs0) =
-   let go constrs curObj ((dir,obj):objs) (Right (Optimal, (opt,_))) =
-         mapFst (opt:) $
-         let extConstrs = (curObj==.opt) : constrs in
-         go extConstrs obj objs $
-         solver extConstrs (dir, objectiveFromTerms sh obj)
-       go _ _ _ sol = ([], sol)
-   in go constrs0 (snd obj0) objs0 $
-      solver constrs0 $ mapSnd (objectiveFromTerms sh) obj0
-
-
-{-# DEPRECATED simplexSuccessive "use GLPK.Monad instead" #-}
-{-# DEPRECATED exactSuccessive "use GLPK.Monad instead" #-}
-{-# DEPRECATED interiorSuccessive "run 'interior' in Either monad instead" #-}
 {- |
-Like the @Multi@ functions,
-but allows not only to fix the previously
-found optimal solution as constraint,
-but allows constraints with a tolerance.
-This is necessary in the presence of rounding errors.
-
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.simplexSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.exactSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False
--}
-simplexSuccessive, exactSuccessive, interiorSuccessive ::
-   (Traversable f, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
-   Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) ->
-   f ((SolutionType, (Double, Array sh Double)) -> Constraints ix,
-      (Direction, Objective sh)) ->
-   Either NoSolutionType
-      (NonEmpty.T f (SolutionType, (Double, Array sh Double)))
-simplexSuccessive = solveSuccessiveInPlace (flip FFI.glp_simplex nullPtr)
-exactSuccessive = solveSuccessiveInPlace (flip FFI.glp_exact nullPtr)
-interiorSuccessive = solveSuccessive . interior
-
-{-# DEPRECATED solveSuccessive
-      "run simple solvers in GLPK.Monad or Either monad instead" #-}
-{- |
-Allows for generic implementation of 'simplexSuccessive' et.al.
-without reuse of interim results.
-
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.simplexSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.simplex bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.exactSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.exact bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB
--}
-solveSuccessive ::
-   (Traversable f, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
-   (Constraints ix -> (Direction, Objective sh) -> Solution sh) ->
-   Constraints ix ->
-   (Direction, Objective sh) ->
-   f ((SolutionType, (Double, Array sh Double)) -> Constraints ix,
-      (Direction, Objective sh)) ->
-   Either NoSolutionType
-      (NonEmpty.T f (SolutionType, (Double, Array sh Double)))
-solveSuccessive solver constrs0 obj0 objs = do
-   let checkShape obj =
-         if Array.shape (snd obj0) == Array.shape obj
-            then obj
-            else error "GLPK.solveSuccessive: objective shapes mismatch"
-   let solveWithConstraints constrs problem =
-         (\sol -> (sol, (constrs,sol))) <$> solver constrs problem
-   (sol0,state0) <- solveWithConstraints constrs0 obj0
-   NonEmpty.cons sol0 <$>
-      MS.evalStateT
-         (for objs $
-            \(newConstrs,(dir,obj)) -> MS.StateT $ \(constrs,sol) ->
-               solveWithConstraints
-                  (newConstrs sol ++ constrs)
-                  (dir, checkShape obj))
-         state0
-
-
-{- |
 >>> case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.exact [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
 Right (Optimal,(28.0,(5.0,0.0,4.0)))
 
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.exact bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.exact bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.exact bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.exact bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.01 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.exact bounds constrs (dir,obj) of Right (LP.Optimal, (opt,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> let val = TestLP.scalarProduct obj $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin) in (case dir of LP.Minimize -> opt-0.01 <= val; LP.Maximize -> opt+0.01 >= val); _ -> QC.property False
 -}
 exact ::
    (Shape.Indexed sh, Shape.Index sh ~ ix) =>
    Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) -> Solution sh
+   (Direction, Objective sh) -> Result sh
 exact = solve (flip FFI.glp_exact nullPtr)
 
 
@@ -285,70 +131,25 @@
    (Shape.Indexed sh, Shape.Index sh ~ ix) =>
    (Foreign.Ptr FFI.Problem -> IO FFI.GlpkSimplexStatus) ->
    Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) -> Solution sh
+   (Direction, Objective sh) -> Result sh
 solve solver bounds constrs (dir,obj) = unsafePerformIO $
    bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> do
    storeProblem bounds constrs (dir,obj) lp
    void $ solver lp
    peekSimplexSolution (Array.shape obj) lp
 
-{-# INLINE solveSuccessiveInPlace #-}
-solveSuccessiveInPlace ::
-   (Traversable f, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
-   (Foreign.Ptr FFI.Problem -> IO FFI.GlpkSimplexStatus) ->
-   Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) ->
-   f ((SolutionType, (Double, Array sh Double)) -> Constraints ix,
-      (Direction, Objective sh)) ->
-   Either NoSolutionType
-      (NonEmpty.T f (SolutionType, (Double, Array sh Double)))
-solveSuccessiveInPlace solver bounds constrs0 (dir0,obj0) objs =
-      unsafePerformIO $
-   bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> ME.runExceptT $ do
-   let shape = Array.shape obj0
-   sol0 <- ME.ExceptT $ do
-      storeProblem bounds constrs0 (dir0,obj0) lp
-      void $ solver lp
-      peekSimplexSolution shape lp
-   NonEmpty.cons sol0 <$>
-      MS.evalStateT
-         (for objs $
-            \(makeNewConstrs,(dir,obj)) -> MS.StateT $ \sol ->
-                  fmap (\sol1 -> (sol1, sol1)) $ ME.ExceptT $ do
-               setDirection lp dir
-               when (shape /= Array.shape obj) $
-                  error "GLPK.solveSuccessiveInplace: objective shapes mismatch"
-               setObjective lp obj
-               let newConstrs = makeNewConstrs sol
-               newRow <- FFI.glp_add_rows lp $ fromIntegral $ length newConstrs
-               for_ (zip [newRow..] (map prepareBounds newConstrs)) $
-                     \(row, (terms,(bnd,lo,up))) -> do
-                  FFI.glp_set_row_bnds lp row bnd lo up
-                  let numTerms = length terms
-                  allocaArray numTerms $ \indicesPtr ->
-                     allocaArray numTerms $ \coeffsPtr -> do
-                     for_ (zip [1..] terms) $
-                        \(k, Term c x) -> do
-                           pokeElem indicesPtr k (columnIndex shape x)
-                           pokeElem coeffsPtr k (realToFrac c)
-                     FFI.glp_set_mat_row lp row
-                        (fromIntegral numTerms) indicesPtr coeffsPtr
-               void $ solver lp
-               peekSimplexSolution shape lp)
-         sol0
 
 
-
 {- |
 >>> case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapPair (round3, Array.toTuple . Array.map round3)) <$> LP.interior [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
 Right (Optimal,(28.0,(5.0,0.0,4.0)))
 
-prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.interior bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.interior bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False
 -}
 interior ::
    (Shape.Indexed sh, Shape.Index sh ~ ix) =>
    Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) -> Solution sh
+   (Direction, Objective sh) -> Result sh
 interior bounds constrs (dir,obj) = unsafePerformIO $
    bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> do
    storeProblem bounds constrs (dir,obj) lp
@@ -369,7 +170,7 @@
    Bounds ix -> Constraints ix ->
    (Direction, Objective sh) -> Foreign.Ptr FFI.Problem -> IO ()
 storeProblem bounds constrs (dir,obj) lp = do
-   void $ FFI.glp_term_out FFI.glpkOff
+   Debug.initLog
    let shape = Array.shape obj
    setDirection lp dir
    firstRow <- FFI.glp_add_rows lp $ fromIntegral $ length constrs
@@ -385,82 +186,8 @@
       for_ (zip [1..] $ concat $
             zipWith (map . (,)) [firstRow..] $
             map (fst . prepareBounds) constrs) $
-         \(k, (row, Term c x)) -> do
+         \(k, (row, LP.Term c x)) -> do
             pokeElem ia k row
             pokeElem ja k (columnIndex shape x)
             pokeElem ar k (realToFrac c)
       FFI.glp_load_matrix lp (fromIntegral numTerms) ia ja ar
-
-
-
-
-class FormatIdentifier ix where
-   formatIdentifier :: ix -> String
-
-instance FormatIdentifier Char where
-   formatIdentifier x = [x]
-
-instance FormatIdentifier c => FormatIdentifier [c] where
-   formatIdentifier = concatMap formatIdentifier
-
-instance FormatIdentifier Int where
-   formatIdentifier = printf "x%d"
-
-instance FormatIdentifier Integer where
-   formatIdentifier = printf "x%d"
-
-
-formatBound :: (FormatIdentifier ix) => Inequality ix -> String
-formatBound (Inequality ix bnd) =
-   printf "var %s%s;" (formatIdentifier ix) $
-   case bnd of
-      LessEqual up -> printf ", <=%f" up
-      GreaterEqual lo -> printf ", >=%f" lo
-      Between lo up -> printf ", >=%f, <=%f" lo up
-      Equal x -> printf ", =%f" x
-      Free -> ""
-
-
-formatSum :: (FormatIdentifier ix) => [Term ix] -> String
-formatSum [] = "0"
-formatSum xs =
-   let formatTerm (Term c ix) = printf "%f*%s" c (formatIdentifier ix) in
-   List.intercalate "+" $ map formatTerm xs
-
-formatConstraint :: (FormatIdentifier ix) => Inequality [Term ix] -> String
-formatConstraint (Inequality terms bnd) =
-   let sumStr = formatSum terms in
-   case bnd of
-      LessEqual up -> printf "%s <= %f" sumStr up
-      GreaterEqual lo -> printf "%f <= %s" lo sumStr
-      Between lo up -> printf "%f <= %s <= %f" lo sumStr up
-      Equal x -> printf "%s = %f" sumStr x
-      Free -> sumStr
-
-formatDirection :: Direction -> String
-formatDirection Minimize = "minimize"
-formatDirection Maximize = "maximize"
-
-formatObjective ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix, FormatIdentifier ix) =>
-   Objective sh -> String
-formatObjective =
-   formatSum . map (\(ix,c) -> Term c ix) . Array.toAssociations
-
-formatMathProg ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix, FormatIdentifier ix) =>
-   Bounds ix -> Constraints ix ->
-   (Direction, Objective sh) -> [String]
-formatMathProg bounds constrs (dir,obj) =
-   map formatBound bounds ++
-   "" :
-   formatDirection dir :
-   printf "value: %s;" (formatObjective obj) :
-   "" :
-   "subject to" :
-   zipWith
-      (\k constr -> printf "constr%d: %s;" k $ formatConstraint constr)
-      [(0::Int)..] constrs ++
-   "" :
-   "end;" :
-   []
diff --git a/src/Numeric/GLPK/Monad.hs b/src/Numeric/GLPK/Monad.hs
--- a/src/Numeric/GLPK/Monad.hs
+++ b/src/Numeric/GLPK/Monad.hs
@@ -4,6 +4,9 @@
 {- |
 The monadic interface to GLPK allows to optimize
 with respect to multiple objectives, successively.
+
+It is not intended to provide a general imperative interface to GLPK
+that manipulates GLPK's state by setters and getters.
 -}
 module Numeric.GLPK.Monad (
    T,
@@ -14,11 +17,13 @@
    ) where
 
 import qualified Math.Programming.Glpk.Header as FFI
+import qualified Numeric.GLPK.Debug as Debug
 import Numeric.GLPK.Private
-         (Term(Term), Constraints, Solution,
+         (Constraints, Result,
           allocaArray, pokeElem, columnIndex, prepareBounds, storeBounds,
           setDirection, setObjective, peekSimplexSolution)
 import Numeric.GLPK (Bounds, Direction(..), Objective)
+import Numeric.LinearProgramming.Common (Term(Term))
 
 import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
@@ -35,14 +40,82 @@
 
 
 {- $setup
+>>> :set -XTypeFamilies
+>>> :set -XTypeOperators
+>>> import qualified Numeric.LinearProgramming.Monad as LPMonad
+>>> import qualified Numeric.LinearProgramming.Test as TestLP
 >>> import qualified Numeric.GLPK.Monad as LP
+>>> import qualified Numeric.GLPK as GLPK
+>>> import Test.Numeric.GLPK.Utility
+>>>    (traverseLag, traverse_Lag, approxSuccession)
 >>> import Test.Numeric.GLPK (TripletShape, tripletShape)
->>> import Numeric.GLPK ((.*), (<=.))
+>>> import Numeric.GLPK (Bounds, Constraints, Objective, (.*), (<=.))
 >>>
 >>> import qualified Data.Array.Comfort.Storable as Array
 >>> import qualified Data.Array.Comfort.Shape as Shape
+>>> import qualified Data.NonEmpty as NonEmpty
+>>> import Data.Array.Comfort.Storable (Array)
+>>> import Data.Traversable (Traversable)
+>>> import Data.Foldable (Foldable)
 >>>
+>>> import qualified Control.Monad.Trans.Except as ME
+>>>
 >>> import Data.Tuple.HT (mapSnd)
+>>>
+>>> import qualified Test.QuickCheck as QC
+>>>
+>>>
+>>> runSuccessive ::
+>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Foldable t) =>
+>>>    sh ->
+>>>    Bounds ix ->
+>>>    (Constraints ix, (GLPK.Direction, Objective sh)) ->
+>>>    t (Double -> Constraints ix, (GLPK.Direction, Objective sh)) ->
+>>>    Either GLPK.FailureType ()
+>>> runSuccessive shape bounds (constrs,dirObj) objs =
+>>>    LP.run shape bounds $ ME.runExceptT $ do
+>>>       (_solType, (opt, _xs)) <- ME.ExceptT $ LP.simplex constrs dirObj
+>>>       traverse_Lag opt
+>>>          (\prevResult (newConstr, dirObjI) -> do
+>>>              (_solType, (optI, _xs)) <-
+>>>                 ME.ExceptT $ LP.simplex (newConstr prevResult) dirObjI
+>>>              return optI)
+>>>          objs
+>>>
+>>> solveSuccessiveWarm ::
+>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>
+>>>    sh ->
+>>>    Bounds ix ->
+>>>    (Constraints ix, (GLPK.Direction, Objective sh)) ->
+>>>    t (Double -> Constraints ix, (GLPK.Direction, Objective sh)) ->
+>>>    Either GLPK.FailureType
+>>>       (NonEmpty.T t (GLPK.SolutionType, (Double, Array sh Double)))
+>>> solveSuccessiveWarm shape bounds (constrs,dirObj) objs =
+>>>    LP.run shape bounds $ ME.runExceptT $ do
+>>>       result <- ME.ExceptT $ LP.simplex constrs dirObj
+>>>       NonEmpty.Cons result <$>
+>>>          traverseLag result
+>>>             (\(_solType, (prevOpt, _xs)) (newConstr, dirObjI) ->
+>>>                 ME.ExceptT $ LP.simplex (newConstr prevOpt) dirObjI)
+>>>             objs
+>>>
+>>> solveSuccessiveGen ::
+>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>
+>>>    sh ->
+>>>    Bounds ix ->
+>>>    (Constraints ix, (GLPK.Direction, Objective sh)) ->
+>>>    t (Double -> Constraints ix, (GLPK.Direction, Objective sh)) ->
+>>>    Either GLPK.FailureType
+>>>       (NonEmpty.T t (GLPK.SolutionType, (Double, Array sh Double)))
+>>> solveSuccessiveGen shape bounds (constrs,dirObj) objs =
+>>>    LPMonad.run shape bounds $ ME.runExceptT $ do
+>>>       result <- ME.ExceptT $ LPMonad.lift GLPK.simplex constrs dirObj
+>>>       NonEmpty.Cons result <$>
+>>>          traverseLag result
+>>>             (\(_solType, (prevOpt, _xs)) (newConstr, dirObjI) ->
+>>>                 ME.ExceptT $
+>>>                    LPMonad.lift GLPK.simplex (newConstr prevOpt) dirObjI)
+>>>             objs
 -}
 
 
@@ -55,7 +128,7 @@
 run shape bounds (Cons act) =
    unsafePerformIO $
    bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> do
-      void $ FFI.glp_term_out FFI.glpkOff
+      Debug.initLog
       storeBounds lp shape bounds
       liftIO $ MR.runReaderT act (shape, lp)
 
@@ -66,17 +139,23 @@
 
 >>> case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd (mapSnd Array.toTuple) <$> LP.run tripletShape [] (LP.simplex [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))
 Right (Optimal,(28.0,(5.0,0.0,4.0)))
+
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (GLPK.simplex bounds constrs (dir,obj), LP.run (Array.shape origin) bounds $ LP.simplex constrs (dir,obj)) of (Right (GLPK.Optimal, (optA,_)), Right (GLPK.Optimal, (optB,_))) -> TestLP.approxReal 0.1 optA optB; _ -> False
+
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> either (\msg -> QC.counterexample (show msg) False) (const $ QC.property True) $ runSuccessive (Array.shape origin) bounds (constrs,dirObj) objs
+
+prop> TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> approxSuccession 0.01 (solveSuccessiveWarm (Array.shape origin) bounds (constrs,dirObj) objs) (solveSuccessiveGen (Array.shape origin) bounds (constrs,dirObj) objs)
 -}
 simplex ::
    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
    Constraints ix ->
-   (Direction, Objective sh) -> T sh (Solution sh)
+   (Direction, Objective sh) -> T sh (Result sh)
 simplex = solve (flip FFI.glp_simplex nullPtr)
 
 exact ::
    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
    Constraints ix ->
-   (Direction, Objective sh) -> T sh (Solution sh)
+   (Direction, Objective sh) -> T sh (Result sh)
 exact = solve (flip FFI.glp_exact nullPtr)
 
 solve ::
@@ -84,7 +163,7 @@
    (Ptr FFI.Problem -> IO FFI.GlpkSimplexStatus) ->
    Constraints ix ->
    (Direction, Objective sh) ->
-   T sh (Solution sh)
+   T sh (Result sh)
 solve method constrs (dir,obj) = Cons $ do
    (shape, lp) <- MR.ask
    when (shape /= Array.shape obj) $
diff --git a/src/Numeric/GLPK/Private.hs b/src/Numeric/GLPK/Private.hs
--- a/src/Numeric/GLPK/Private.hs
+++ b/src/Numeric/GLPK/Private.hs
@@ -2,6 +2,11 @@
 {-# LANGUAGE TypeOperators #-}
 module Numeric.GLPK.Private where
 
+import qualified Numeric.LinearProgramming.Common as LP
+import Numeric.LinearProgramming.Common
+         (Bound(..), Inequality(Inequality),
+          Bounds, Direction(..), Objective)
+
 import qualified Math.Programming.Glpk.Header as FFI
 
 import qualified Data.Array.Comfort.Storable.Mutable as Mutable
@@ -18,36 +23,11 @@
 import Foreign.C.Types (CDouble)
 
 
-
-data Term ix = Term Double ix
-   deriving (Show)
-
-
-data Inequality x = Inequality x Bound
-   deriving Show
-
-data Bound =
-     LessEqual Double
-   | GreaterEqual Double
-   | Between Double Double
-   | Equal Double
-   | Free
-   deriving Show
-
-instance Functor Inequality where
-   fmap f (Inequality x bnd)  =  Inequality (f x) bnd
-
-type Bounds ix = [Inequality ix]
-
-type Constraints ix = [Inequality [Term ix]]
-
-data Direction = Minimize | Maximize
-   deriving (Eq, Enum, Bounded, Show)
-
-type Objective sh = Array sh Double
+type Term = LP.Term Double
 
+type Constraints ix = LP.Constraints Double ix
 
-data NoSolutionType =
+data FailureType =
      Undefined
    | NoFeasible
    | Unbounded
@@ -59,7 +39,7 @@
    | Optimal
    deriving (Eq, Show)
 
-instance NFData NoSolutionType where
+instance NFData FailureType where
     rnf NoFeasible = ()
     rnf _ = ()
 
@@ -67,8 +47,8 @@
     rnf Optimal = ()
     rnf _ = ()
 
-type Solution sh =
-      Either NoSolutionType (SolutionType, (Double, Array sh Double))
+type Result sh =
+      Either FailureType (SolutionType, (Double, Array sh Double))
 
 
 {- |
@@ -151,7 +131,7 @@
    for_ (Shape.indices defShape) (act arr)
    Array.reshape shape <$> Mutable.freeze arr
 
-analyzeStatus :: FFI.GlpkSolutionStatus -> Either NoSolutionType SolutionType
+analyzeStatus :: FFI.GlpkSolutionStatus -> Either FailureType SolutionType
 analyzeStatus status =
    fromMaybe (error "glpk-solver: solution type unknown") $ lookup status $
       (FFI.glpkUndefined,  Left Undefined) :
@@ -164,7 +144,7 @@
 
 
 peekSimplexSolution ::
-   (Shape.C sh) => sh -> Foreign.Ptr FFI.Problem -> IO (Solution sh)
+   (Shape.C sh) => sh -> Foreign.Ptr FFI.Problem -> IO (Result sh)
 peekSimplexSolution shape lp = do
    let examine =
          liftA2 (,)
diff --git a/src/debug-off/Numeric/GLPK/Debug.hs b/src/debug-off/Numeric/GLPK/Debug.hs
new file mode 100644
--- /dev/null
+++ b/src/debug-off/Numeric/GLPK/Debug.hs
@@ -0,0 +1,8 @@
+module Numeric.GLPK.Debug where
+
+import qualified Math.Programming.Glpk.Header as FFI
+import Control.Monad (void)
+
+
+initLog :: IO ()
+initLog = void $ FFI.glp_term_out FFI.glpkOff
diff --git a/src/debug-on/Numeric/GLPK/Debug.hs b/src/debug-on/Numeric/GLPK/Debug.hs
new file mode 100644
--- /dev/null
+++ b/src/debug-on/Numeric/GLPK/Debug.hs
@@ -0,0 +1,5 @@
+module Numeric.GLPK.Debug where
+
+
+initLog :: IO ()
+initLog = return ()
diff --git a/test/Test/Numeric/GLPK.hs b/test/Test/Numeric/GLPK.hs
--- a/test/Test/Numeric/GLPK.hs
+++ b/test/Test/Numeric/GLPK.hs
@@ -1,13 +1,13 @@
 -- Do not edit! Automatically created with doctest-extract from src/Numeric/GLPK.hs
-{-# LINE 82 "src/Numeric/GLPK.hs" #-}
+{-# LINE 67 "src/Numeric/GLPK.hs" #-}
 
 module Test.Numeric.GLPK where
 
 import Test.DocTest.Base
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 83 "src/Numeric/GLPK.hs" #-}
-import     qualified Test.Numeric.GLPK.Generator as TestLP
+{-# LINE 68 "src/Numeric/GLPK.hs" #-}
+import     qualified Numeric.LinearProgramming.Test as TestLP
 import     qualified Numeric.GLPK as LP
 import     Numeric.GLPK ((.*), (<=.), (==.))
 
@@ -34,71 +34,76 @@
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Numeric.GLPK:141: "
-{-# LINE 141 "src/Numeric/GLPK.hs" #-}
+ DocTest.printPrefix "Numeric.GLPK:99: "
+{-# LINE 99 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 141 "src/Numeric/GLPK.hs" #-}
+{-# LINE 99 "src/Numeric/GLPK.hs" #-}
      (\target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (1,1) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (TestLP.approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False)
- DocTest.printPrefix "Numeric.GLPK:142: "
-{-# LINE 142 "src/Numeric/GLPK.hs" #-}
+ DocTest.printPrefix "Numeric.GLPK:100: "
+{-# LINE 100 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 142 "src/Numeric/GLPK.hs" #-}
+{-# LINE 100 "src/Numeric/GLPK.hs" #-}
      (\(QC.Positive posWeight) (QC.Positive negWeight) target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (posWeight,negWeight) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False)
- DocTest.printPrefix "Numeric.GLPK:143: "
-{-# LINE 143 "src/Numeric/GLPK.hs" #-}
+ DocTest.printPrefix "Numeric.GLPK:101: "
+{-# LINE 101 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 143 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, _) -> True; _ -> False)
- DocTest.printPrefix "Numeric.GLPK:138: "
-{-# LINE 138 "src/Numeric/GLPK.hs" #-}
+{-# LINE 101 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, _) -> True; _ -> False)
+ DocTest.printPrefix "Numeric.GLPK:102: "
+{-# LINE 102 "src/Numeric/GLPK.hs" #-}
+ DocTest.property
+{-# LINE 102 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:103: "
+{-# LINE 103 "src/Numeric/GLPK.hs" #-}
+ DocTest.property
+{-# LINE 103 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.1 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:104: "
+{-# LINE 104 "src/Numeric/GLPK.hs" #-}
+ DocTest.property
+{-# LINE 104 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, (opt,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> let val = TestLP.scalarProduct obj $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin) in (case dir of LP.Minimize -> opt-0.01 <= val; LP.Maximize -> opt+0.01 >= val); _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:96: "
+{-# LINE 96 "src/Numeric/GLPK.hs" #-}
  DocTest.example
-{-# LINE 138 "src/Numeric/GLPK.hs" #-}
+{-# LINE 96 "src/Numeric/GLPK.hs" #-}
    (case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))
   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]
- DocTest.printPrefix "Numeric.GLPK:180: "
-{-# LINE 180 "src/Numeric/GLPK.hs" #-}
- DocTest.property
-{-# LINE 180 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case LP.simplexMulti bounds constrs (Array.shape origin) objs of (_, Right (LP.Optimal, _)) -> QC.property True; result -> QC.counterexample (show result) False)
- DocTest.printPrefix "Numeric.GLPK:219: "
-{-# LINE 219 "src/Numeric/GLPK.hs" #-}
- DocTest.property
-{-# LINE 219 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.simplexSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False)
- DocTest.printPrefix "Numeric.GLPK:220: "
-{-# LINE 220 "src/Numeric/GLPK.hs" #-}
+ DocTest.printPrefix "Numeric.GLPK:117: "
+{-# LINE 117 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 220 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.exactSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False)
- DocTest.printPrefix "Numeric.GLPK:240: "
-{-# LINE 240 "src/Numeric/GLPK.hs" #-}
+{-# LINE 117 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.exact bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:118: "
+{-# LINE 118 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 240 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.simplexSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.simplex bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB)
- DocTest.printPrefix "Numeric.GLPK:241: "
-{-# LINE 241 "src/Numeric/GLPK.hs" #-}
+{-# LINE 118 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.exact bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:119: "
+{-# LINE 119 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 241 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.exactSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.exact bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB)
- DocTest.printPrefix "Numeric.GLPK:274: "
-{-# LINE 274 "src/Numeric/GLPK.hs" #-}
+{-# LINE 119 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.exact bounds constrs (dir,obj) of Right (LP.Optimal, (_,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.01 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:120: "
+{-# LINE 120 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 274 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.exact bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False)
- DocTest.printPrefix "Numeric.GLPK:271: "
-{-# LINE 271 "src/Numeric/GLPK.hs" #-}
+{-# LINE 120 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.exact bounds constrs (dir,obj) of Right (LP.Optimal, (opt,sol)) -> QC.forAll (QC.choose (0,1)) $ \lambda -> let val = TestLP.scalarProduct obj $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin) in (case dir of LP.Minimize -> opt-0.01 <= val; LP.Maximize -> opt+0.01 >= val); _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:114: "
+{-# LINE 114 "src/Numeric/GLPK.hs" #-}
  DocTest.example
-{-# LINE 271 "src/Numeric/GLPK.hs" #-}
+{-# LINE 114 "src/Numeric/GLPK.hs" #-}
    (case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.exact [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))
   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]
- DocTest.printPrefix "Numeric.GLPK:346: "
-{-# LINE 346 "src/Numeric/GLPK.hs" #-}
+ DocTest.printPrefix "Numeric.GLPK:147: "
+{-# LINE 147 "src/Numeric/GLPK.hs" #-}
  DocTest.property
-{-# LINE 346 "src/Numeric/GLPK.hs" #-}
-     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.interior bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False)
- DocTest.printPrefix "Numeric.GLPK:343: "
-{-# LINE 343 "src/Numeric/GLPK.hs" #-}
+{-# LINE 147 "src/Numeric/GLPK.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.interior bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False)
+ DocTest.printPrefix "Numeric.GLPK:144: "
+{-# LINE 144 "src/Numeric/GLPK.hs" #-}
  DocTest.example
-{-# LINE 343 "src/Numeric/GLPK.hs" #-}
+{-# LINE 144 "src/Numeric/GLPK.hs" #-}
    (case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapPair (round3, Array.toTuple . Array.map round3)) <$> LP.interior [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))
   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]
diff --git a/test/Test/Numeric/GLPK/Generator.hs b/test/Test/Numeric/GLPK/Generator.hs
deleted file mode 100644
--- a/test/Test/Numeric/GLPK/Generator.hs
+++ /dev/null
@@ -1,208 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-module Test.Numeric.GLPK.Generator where
-
-import qualified Numeric.GLPK as LP
-import Numeric.GLPK ((<=.), (>=.))
-
-import qualified Test.QuickCheck as QC
-import System.Random (Random)
-
-import qualified Data.Array.Comfort.Boxed as BoxedArray
-import qualified Data.Array.Comfort.Storable as Array
-import qualified Data.Array.Comfort.Shape as Shape
-import qualified Data.NonEmpty as NonEmpty
-import qualified Data.List.HT as ListHT
-import qualified Data.Ix as Ix
-import Data.Array.Comfort.Storable (Array, (!))
-import Data.Traversable (sequenceA, for)
-import Data.Tuple.HT (mapSnd)
-import Data.Maybe (fromMaybe)
-import Data.Int (Int64)
-
-import Control.Applicative (liftA2)
-
-import Text.Printf (PrintfArg, printf)
-
-import Foreign.Storable (Storable)
-
-
-{- |
-Generate constraints in the form of a polyhedron
-which contains warrantedly the zero vector.
-That is, there is an admissible solution.
-In order to assert that the polyhedron is closed,
-we bound all variables by a hypercube.
--}
-genProblem ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
-   Array sh a -> QC.Gen (LP.Bounds ix, LP.Constraints ix)
-genProblem origin =
-   liftA2 (,)
-      (for (Array.toAssociations origin) $ \(ix,x) ->
-         LP.Inequality ix <$>
-         liftA2 LP.Between
-            (doubleFromElement . (x+) <$> QC.choose (-100,-50))
-            (doubleFromElement . (x+) <$> QC.choose (50,100)))
-      (do
-         numConstraints <- QC.choose (1,20)
-         QC.vectorOf numConstraints $ do
-            ixs <- QC.sublistOf $ Shape.indices $ Array.shape origin
-            terms <- for ixs $ \ix -> do
-               coeff <- QC.choose (-10,10)
-               return (coeff, ix)
-            let offset = scalarProduct terms origin
-            let deviation = 25
-            LP.Inequality
-               (map (uncurry (LP.Term . doubleFromElement)) terms)
-               <$>
-               QC.oneof (
-                  (do bound <- QC.choose (offset-deviation, offset+deviation)
-                      return $
-                         if bound > offset
-                            then LP.LessEqual    $ doubleFromElement bound
-                            else LP.GreaterEqual $ doubleFromElement bound) :
-                  (liftA2 LP.Between
-                     (doubleFromElement <$>
-                        QC.choose (offset-deviation, offset))
-                     (doubleFromElement <$>
-                        QC.choose (offset, offset+deviation))) :
-                  []))
-
-scalarProduct ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix, Storable a, Num a) =>
-   [(a,ix)] -> Array sh a -> a
-scalarProduct terms origin =
-   sum $ map (\(coeff, ix) -> coeff * origin!ix) terms
-
-genVarShape :: QC.Gen (Shape.Range Char)
-genVarShape = Shape.Range 'a' <$> QC.choose ('a','j')
-
-genOrigin :: QC.Gen (Array (Shape.Range Char) Int64)
-genOrigin = genVector =<< genVarShape
-
-_genOrigin :: QC.Gen (Array (Shape.Range Char) Double)
-_genOrigin = genVector =<< genVarShape
-
-
-shrinkVarShape :: Shape.Range Char -> [Shape.Range Char]
-shrinkVarShape (Shape.Range from to) =
-   if from<to then [Shape.Range from (pred to)] else []
-
-shrinkOrigin ::
-   (Storable a) => Array (Shape.Range Char) a -> [Array (Shape.Range Char) a]
-shrinkOrigin vec =
-   case Array.shape vec of
-      Shape.Range from to ->
-         if from<to
-            then [Array.sample (Shape.Range from (pred to)) (vec!)]
-            else []
-
-
-class (Storable a, Random a, Num a, Ord a) => Element a where
-   doubleFromElement :: a -> Double
-
-instance Element Double where
-   doubleFromElement = id
-
-instance Element Int64 where
-   doubleFromElement = fromIntegral
-
-genObjective ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
-   Array sh a -> QC.Gen (LP.Direction, LP.Objective sh)
-genObjective origin =
-   liftA2 (,) QC.arbitraryBoundedEnum
-      (fmap (Array.map doubleFromElement . flip asTypeOf origin) $
-       genVector $ Array.shape origin)
-
-genVector :: (Shape.Indexed sh, Element a) => sh -> QC.Gen (Array sh a)
-genVector shape =
-   fmap Array.fromBoxed $ sequenceA $
-   BoxedArray.fromAssociations (QC.choose (-10,10)) shape []
-
-shrinkProblem ::
-   (LP.Bounds ix, LP.Constraints ix) ->
-   [(LP.Bounds ix, LP.Constraints ix)]
-shrinkProblem (bounds, constraints) =
-   map (\shrinked -> (bounds, shrinked)) $
-   filter (not . null) $ QC.shrinkList (const []) constraints
-
-genObjectives ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
-   Array sh a -> QC.Gen (NonEmpty.T [] (LP.Direction, [LP.Term ix]))
-genObjectives origin = do
-   let shape = Array.shape origin
-   let stageRange :: (Int,Int)
-       stageRange = (0,3)
-   stages <- for (Shape.indices shape) $ \ix -> (,) ix <$> QC.choose stageRange
-   let varSets =
-         fromMaybe (error "there should be at least one stage") $
-         NonEmpty.fetch $
-         filter (not . null) $
-         map (\k -> map fst $ filter ((k==) . snd) stages) $
-         Ix.range stageRange
-   let asTypeOfElement :: a -> f a -> a
-       asTypeOfElement = const
-   for varSets $ \varSet ->
-      liftA2 (,)
-         QC.arbitraryBoundedEnum
-         (for varSet $ \ix ->
-            flip LP.Term ix . doubleFromElement
-               <$> QC.choose (-10, 10 `asTypeOfElement` origin))
-
-shrinkObjectives ::
-   NonEmpty.T [] (LP.Direction, [LP.Term ix]) ->
-   [NonEmpty.T [] (LP.Direction, [LP.Term ix])]
-shrinkObjectives (NonEmpty.Cons obj objs) =
-   map (NonEmpty.Cons obj) $
-   QC.shrinkList
-      (\(dir,terms) ->
-         map ((,) dir) $ filter (not . null) $
-         QC.shrinkList (const []) terms)
-      objs
-
-successiveObjectives ::
-   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
-   Array sh a -> Double ->
-   NonEmpty.T [] (LP.Direction, [LP.Term ix]) ->
-   ((LP.Direction, LP.Objective sh),
-    [((LP.SolutionType, (Double, Array sh Double)) -> LP.Constraints ix,
-      (LP.Direction, LP.Objective sh))])
-successiveObjectives origin tol xs =
-   let shape = Array.shape origin in
-   (mapSnd (LP.objectiveFromTerms shape) $ NonEmpty.head xs,
-    NonEmpty.mapAdjacent
-      (\(dir0,obj0) y1 ->
-         (\(_sol,(opt,_vec)) ->
-            case dir0 of
-               LP.Minimize -> [obj0 <=. opt + tol]
-               LP.Maximize -> [obj0 >=. opt - tol],
-          mapSnd (LP.objectiveFromTerms shape) y1))
-      xs)
-
-
-approxReal :: (Ord a, Num a) => a -> a -> a -> Bool
-approxReal tol x y = abs (x-y) <= tol
-
-approx :: (PrintfArg a, Ord a, Num a) => String -> a -> a -> a -> QC.Property
-approx name tol x y =
-   QC.counterexample (printf "%s: %f - %f" name x y) (approxReal tol x y)
-
-approxSuccession ::
-   (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>
-   a ->
-   Either LP.NoSolutionType
-      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->
-   Either LP.NoSolutionType
-      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->
-   QC.Property
-approxSuccession tol x y =
-   QC.counterexample (show x) $
-   QC.counterexample (show y) $
-   case (x,y) of
-      (Left sx, Left sy) -> sx==sy
-      (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->
-         let equalSol (solX, (optX, _)) (solY, (optY, _)) =
-               solX == solY && approxReal tol optX optY
-         in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys
-      _ -> False
diff --git a/test/Test/Numeric/GLPK/Monad.hs b/test/Test/Numeric/GLPK/Monad.hs
--- a/test/Test/Numeric/GLPK/Monad.hs
+++ b/test/Test/Numeric/GLPK/Monad.hs
@@ -1,26 +1,109 @@
 -- Do not edit! Automatically created with doctest-extract from src/Numeric/GLPK/Monad.hs
-{-# LINE 37 "src/Numeric/GLPK/Monad.hs" #-}
+{-# LINE 42 "src/Numeric/GLPK/Monad.hs" #-}
 
+{-# OPTIONS_GHC -XTypeFamilies #-}
+{-# OPTIONS_GHC -XTypeOperators #-}
 module Test.Numeric.GLPK.Monad where
 
 import Test.DocTest.Base
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 38 "src/Numeric/GLPK/Monad.hs" #-}
+{-# LINE 45 "src/Numeric/GLPK/Monad.hs" #-}
+import     qualified Numeric.LinearProgramming.Monad as LPMonad
+import     qualified Numeric.LinearProgramming.Test as TestLP
 import     qualified Numeric.GLPK.Monad as LP
+import     qualified Numeric.GLPK as GLPK
+import     Test.Numeric.GLPK.Utility
+       (traverseLag, traverse_Lag, approxSuccession)
 import     Test.Numeric.GLPK (TripletShape, tripletShape)
-import     Numeric.GLPK ((.*), (<=.))
+import     Numeric.GLPK (Bounds, Constraints, Objective, (.*), (<=.))
 
 import     qualified Data.Array.Comfort.Storable as Array
 import     qualified Data.Array.Comfort.Shape as Shape
+import     qualified Data.NonEmpty as NonEmpty
+import     Data.Array.Comfort.Storable (Array)
+import     Data.Traversable (Traversable)
+import     Data.Foldable (Foldable)
 
+import     qualified Control.Monad.Trans.Except as ME
+
 import     Data.Tuple.HT (mapSnd)
 
+import     qualified Test.QuickCheck as QC
+
+
+runSuccessive     ::
+       (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Foldable t) =>
+       sh ->
+       Bounds ix ->
+       (Constraints ix, (GLPK.Direction, Objective sh)) ->
+       t (Double -> Constraints ix, (GLPK.Direction, Objective sh)) ->
+       Either GLPK.FailureType ()
+runSuccessive     shape bounds (constrs,dirObj) objs =
+       LP.run shape bounds $ ME.runExceptT $ do
+          (_solType, (opt, _xs)) <- ME.ExceptT $ LP.simplex constrs dirObj
+          traverse_Lag opt
+             (\prevResult (newConstr, dirObjI) -> do
+                 (_solType, (optI, _xs)) <-
+                    ME.ExceptT $ LP.simplex (newConstr prevResult) dirObjI
+                 return optI)
+             objs
+
+solveSuccessiveWarm     ::
+       (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>
+       sh ->
+       Bounds ix ->
+       (Constraints ix, (GLPK.Direction, Objective sh)) ->
+       t (Double -> Constraints ix, (GLPK.Direction, Objective sh)) ->
+       Either GLPK.FailureType
+          (NonEmpty.T t (GLPK.SolutionType, (Double, Array sh Double)))
+solveSuccessiveWarm     shape bounds (constrs,dirObj) objs =
+       LP.run shape bounds $ ME.runExceptT $ do
+          result <- ME.ExceptT $ LP.simplex constrs dirObj
+          NonEmpty.Cons result <$>
+             traverseLag result
+                (\(_solType, (prevOpt, _xs)) (newConstr, dirObjI) ->
+                    ME.ExceptT $ LP.simplex (newConstr prevOpt) dirObjI)
+                objs
+
+solveSuccessiveGen     ::
+       (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>
+       sh ->
+       Bounds ix ->
+       (Constraints ix, (GLPK.Direction, Objective sh)) ->
+       t (Double -> Constraints ix, (GLPK.Direction, Objective sh)) ->
+       Either GLPK.FailureType
+          (NonEmpty.T t (GLPK.SolutionType, (Double, Array sh Double)))
+solveSuccessiveGen     shape bounds (constrs,dirObj) objs =
+       LPMonad.run shape bounds $ ME.runExceptT $ do
+          result <- ME.ExceptT $ LPMonad.lift GLPK.simplex constrs dirObj
+          NonEmpty.Cons result <$>
+             traverseLag result
+                (\(_solType, (prevOpt, _xs)) (newConstr, dirObjI) ->
+                    ME.ExceptT $
+                       LPMonad.lift GLPK.simplex (newConstr prevOpt) dirObjI)
+                objs
+
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Numeric.GLPK.Monad:67: "
-{-# LINE 67 "src/Numeric/GLPK/Monad.hs" #-}
+ DocTest.printPrefix "Numeric.GLPK.Monad:143: "
+{-# LINE 143 "src/Numeric/GLPK/Monad.hs" #-}
+ DocTest.property
+{-# LINE 143 "src/Numeric/GLPK/Monad.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (GLPK.simplex bounds constrs (dir,obj), LP.run (Array.shape origin) bounds $ LP.simplex constrs (dir,obj)) of (Right (GLPK.Optimal, (optA,_)), Right (GLPK.Optimal, (optB,_))) -> TestLP.approxReal 0.1 optA optB; _ -> False)
+ DocTest.printPrefix "Numeric.GLPK.Monad:145: "
+{-# LINE 145 "src/Numeric/GLPK/Monad.hs" #-}
+ DocTest.property
+{-# LINE 145 "src/Numeric/GLPK/Monad.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> either (\msg -> QC.counterexample (show msg) False) (const $ QC.property True) $ runSuccessive (Array.shape origin) bounds (constrs,dirObj) objs)
+ DocTest.printPrefix "Numeric.GLPK.Monad:147: "
+{-# LINE 147 "src/Numeric/GLPK/Monad.hs" #-}
+ DocTest.property
+{-# LINE 147 "src/Numeric/GLPK/Monad.hs" #-}
+     (TestLP.forAllOrigin $ \origin -> TestLP.forAllProblem origin $ \bounds constrs -> TestLP.forAllObjectives origin $ \objs_ -> case TestLP.successiveObjectives origin 0.01 objs_ of (dirObj, objs) -> approxSuccession 0.01 (solveSuccessiveWarm (Array.shape origin) bounds (constrs,dirObj) objs) (solveSuccessiveGen (Array.shape origin) bounds (constrs,dirObj) objs))
+ DocTest.printPrefix "Numeric.GLPK.Monad:140: "
+{-# LINE 140 "src/Numeric/GLPK/Monad.hs" #-}
  DocTest.example
-{-# LINE 67 "src/Numeric/GLPK/Monad.hs" #-}
+{-# LINE 140 "src/Numeric/GLPK/Monad.hs" #-}
    (case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd (mapSnd Array.toTuple) <$> LP.run tripletShape [] (LP.simplex [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)))
   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]
diff --git a/test/Test/Numeric/GLPK/Utility.hs b/test/Test/Numeric/GLPK/Utility.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/GLPK/Utility.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Test.Numeric.GLPK.Utility where
+
+import qualified Numeric.GLPK as LP
+import Numeric.LinearProgramming.Test (approxReal)
+
+import qualified Test.QuickCheck as QC
+
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.NonEmpty as NonEmpty
+import qualified Data.List.HT as ListHT
+import Data.Array.Comfort.Storable (Array)
+import Data.Tuple.HT (double)
+import Data.Traversable (Traversable, traverse)
+import Data.Foldable (Foldable, traverse_)
+
+import qualified Control.Monad.Trans.State as MS
+
+import Foreign.Storable (Storable)
+
+
+approxSuccession ::
+   (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>
+   a ->
+   Either LP.FailureType
+      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->
+   Either LP.FailureType
+      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->
+   QC.Property
+approxSuccession tol x y =
+   QC.counterexample (show x) $
+   QC.counterexample (show y) $
+   case (x,y) of
+      (Left sx, Left sy) -> sx==sy
+      (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->
+         let equalSol (solX, (optX, _)) (solY, (optY, _)) =
+               solX == solY && approxReal tol optX optY
+         in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys
+      _ -> False
+
+
+traverse_Lag ::
+   (Foldable t, Monad m) =>
+   b -> (b -> a -> m b) -> t a -> m ()
+traverse_Lag b0 f =
+   flip MS.evalStateT b0 .
+   traverse_ (\a -> MS.StateT $ \b -> fmap double $ f b a)
+
+traverseLag ::
+   (Traversable t, Monad m) =>
+   b -> (b -> a -> m b) -> t a -> m (t b)
+traverseLag b0 f =
+   flip MS.evalStateT b0 .
+   traverse (\a -> MS.StateT $ \b -> fmap double $ f b a)