diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Henning Thielemann 2023
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Makefile b/Makefile
new file mode 100644
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,4 @@
+run-test:
+	runhaskell Setup configure --user --enable-tests
+	runhaskell Setup build
+	runhaskell Setup haddock
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/linear-programming.cabal b/linear-programming.cabal
new file mode 100644
--- /dev/null
+++ b/linear-programming.cabal
@@ -0,0 +1,44 @@
+Cabal-Version:    2.2
+Name:             linear-programming
+Version:          0.0
+License:          BSD-3-Clause
+License-File:     LICENSE
+Author:           Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Category:         Math
+Tested-With:      GHC ==8.6.5
+Build-Type:       Simple
+Synopsis:         Linear Programming basic definitions
+Description:
+  Basic types and generic functions for use in the packages
+  @coinor-clp@ and @comfort-glpk@.
+Extra-Source-Files:
+  Makefile
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    https://hub.darcs.net/thielema/linear-programming/
+
+Source-Repository head
+  Type:        darcs
+  Location:    https://hub.darcs.net/thielema/linear-programming/
+
+Library
+  Build-Depends:
+    comfort-array >=0.4 && <0.6,
+    QuickCheck >=2.1 && <3,
+    random >=1.0 && <1.3,
+    transformers >=0.3 && <0.7,
+    non-empty >=0.3.2 && <0.4,
+    utility-ht >=0.0.16 && <0.1,
+    base >=4.5 && <5
+
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   src
+  Default-Language: Haskell98
+  Exposed-Modules:
+    Numeric.LinearProgramming.Common
+    Numeric.LinearProgramming.Format
+    Numeric.LinearProgramming.Monad
+    Numeric.LinearProgramming.Test
diff --git a/src/Numeric/LinearProgramming/Common.hs b/src/Numeric/LinearProgramming/Common.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearProgramming/Common.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Numeric.LinearProgramming.Common (
+   Term(..), (.*),
+   Inequality(..),
+   Bound(..),
+   Bounds,
+   Constraints,
+   Direction(..),
+   Objective,
+   free, (<=.), (>=.), (==.), (>=<.),
+   objectiveFromTerms,
+   ) where
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable (Array)
+
+
+
+data Term a ix = Term a ix
+   deriving (Show)
+
+
+infix 7 .*
+
+(.*) :: a -> ix -> Term a ix
+(.*) = Term
+
+
+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 a ix = [Inequality [Term a ix]]
+
+data Direction = Minimize | Maximize
+   deriving (Eq, Enum, Bounded, Show)
+
+type Objective sh = Array sh Double
+
+
+
+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 Double ix] -> Objective sh
+objectiveFromTerms sh =
+   Array.fromAssociations 0 sh . map (\(Term x ix) -> (ix,x))
diff --git a/src/Numeric/LinearProgramming/Format.hs b/src/Numeric/LinearProgramming/Format.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearProgramming/Format.hs
@@ -0,0 +1,97 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Numeric.LinearProgramming.Format (
+   Identifier,
+   mathProg,
+   ) where
+
+import qualified Numeric.LinearProgramming.Common as LP
+import Numeric.LinearProgramming.Common
+         (Bound(..), Inequality(Inequality),
+          Bounds, Direction(..), Objective, (.*))
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.List as List
+
+import Text.Printf (printf)
+
+import Prelude hiding (sum)
+
+
+
+type Term = LP.Term Double
+
+type Constraints ix = LP.Constraints Double ix
+
+
+class Identifier ix where
+   identifier :: ix -> String
+
+instance Identifier Char where
+   identifier x = [x]
+
+instance Identifier c => Identifier [c] where
+   identifier = concatMap identifier
+
+instance Identifier Int where
+   identifier = printf "x%d"
+
+instance Identifier Integer where
+   identifier = printf "x%d"
+
+
+bound :: (Identifier ix) => Inequality ix -> String
+bound (Inequality ix bnd) =
+   printf "var %s%s;" (identifier 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 -> ""
+
+
+sum :: (Identifier ix) => [Term ix] -> String
+sum [] = "0"
+sum xs =
+   let formatTerm (LP.Term c ix) = printf "%f*%s" c (identifier ix) in
+   List.intercalate "+" $ map formatTerm xs
+
+constraint :: (Identifier ix) => Inequality [Term ix] -> String
+constraint (Inequality terms bnd) =
+   let sumStr = sum 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
+
+direction :: Direction -> String
+direction Minimize = "minimize"
+direction Maximize = "maximize"
+
+objective ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Identifier ix) =>
+   Objective sh -> String
+objective =
+   sum . map (\(ix,c) -> c .* ix) . Array.toAssociations
+
+mathProg ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Identifier ix) =>
+   Bounds ix -> Constraints ix ->
+   (Direction, Objective sh) -> [String]
+mathProg bounds constrs (dir,obj) =
+   map bound bounds ++
+   "" :
+   direction dir :
+   printf "value: %s;" (objective obj) :
+   "" :
+   "subject to" :
+   zipWith
+      (\k constr -> printf "constr%d: %s;" k $ constraint constr)
+      [(0::Int)..] constrs ++
+   "" :
+   "end;" :
+   []
diff --git a/src/Numeric/LinearProgramming/Monad.hs b/src/Numeric/LinearProgramming/Monad.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearProgramming/Monad.hs
@@ -0,0 +1,51 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{- |
+Generic implementation of a monad that collects constraints
+over multiple stages.
+It can be used to test solvers that allow for warm start
+or for solvers that do not allow for warm start at all
+(like GLPK's interior point solver).
+-}
+module Numeric.LinearProgramming.Monad (
+   T,
+   run,
+   lift,
+   ) where
+
+import Numeric.LinearProgramming.Common
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+
+import qualified Control.Monad.Trans.RWS as MRWS
+import Control.Monad (when)
+
+
+newtype T sh a =
+   Cons (MRWS.RWS
+            (sh, Bounds (Shape.Index sh))
+            ()
+            (Constraints Double (Shape.Index sh))
+            a)
+      deriving (Functor, Applicative, Monad)
+
+
+run ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
+   sh -> Bounds ix -> T sh a -> a
+run shape bounds (Cons act) =
+   fst $ MRWS.evalRWS act (shape, bounds) []
+
+lift ::
+   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
+   (Bounds ix -> Constraints Double ix -> (Direction, Objective sh) -> a) ->
+   Constraints Double ix -> (Direction, Objective sh) -> T sh a
+lift solver constrs dirObj@(_dir,obj) = Cons $ do
+   (shape,bounds) <- MRWS.ask
+   when (shape /= Array.shape obj) $
+      error "LinearProgramming.Monad.solve: objective shape mismatch"
+   MRWS.modify (constrs++)
+   allConstrs <- MRWS.get
+   return $ solver bounds allConstrs dirObj
diff --git a/src/Numeric/LinearProgramming/Test.hs b/src/Numeric/LinearProgramming/Test.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LinearProgramming/Test.hs
@@ -0,0 +1,282 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Numeric.LinearProgramming.Test (
+   Element,
+   forAllOrigin,
+   forAllProblem,
+   genObjective,
+   forAllObjectives,
+   successiveObjectives,
+   approxReal,
+   approx,
+   checkFeasibility,
+   affineCombination,
+   scalarProduct,
+   ) where
+
+import qualified Numeric.LinearProgramming.Common as LP
+import Numeric.LinearProgramming.Common ((<=.), (>=.), (.*))
+
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck ((.&&.))
+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.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)
+
+
+
+type Term = LP.Term Double
+type Constraints ix = LP.Constraints Double ix
+
+
+{- |
+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, 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 = scalarProductTerms terms origin
+            let deviation = 25
+            LP.Inequality
+               (map (uncurry ((.*) . 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))) :
+                  []))
+
+scalarProductTerms ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Storable a, Num a) =>
+   [(a,ix)] -> Array sh a -> a
+scalarProductTerms 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 []
+
+
+forAllOrigin ::
+   (QC.Testable prop) =>
+   (Array (Shape.Range Char) Int64 -> prop) -> QC.Property
+forAllOrigin = QC.forAllShrink genOrigin shrinkOrigin
+
+
+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 []
+--    BoxedArray.constant shape (QC.choose (-10,10))
+
+shrinkProblem ::
+   (LP.Bounds ix, Constraints ix) ->
+   [(LP.Bounds ix, Constraints ix)]
+shrinkProblem (bounds, constraints) =
+   map (\shrinked -> (bounds, shrinked)) $
+   filter (not . null) $ QC.shrinkList (const []) constraints
+
+forAllProblem ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>
+   (QC.Testable prop, Element a) =>
+   Array sh a -> (LP.Bounds ix -> Constraints ix -> prop) -> QC.Property
+forAllProblem origin =
+   QC.forAllShrink (genProblem origin) shrinkProblem . uncurry
+
+
+genObjectives ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
+   Array sh a -> QC.Gen (NonEmpty.T [] (LP.Direction, [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 ->
+            (.*ix) . doubleFromElement
+               <$> QC.choose (-10, 10 `asTypeOfElement` origin))
+
+shrinkObjectives ::
+   NonEmpty.T [] (LP.Direction, [Term ix]) ->
+   [NonEmpty.T [] (LP.Direction, [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
+
+forAllObjectives ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix, Show ix) =>
+   (QC.Testable prop, Element a) =>
+   Array sh a ->
+   (NonEmpty.T [] (LP.Direction, [Term (Shape.Index sh)]) -> prop) ->
+   QC.Property
+forAllObjectives origin =
+   QC.forAllShrink (genObjectives origin) shrinkObjectives
+
+constraintsFromSolution ::
+   Double -> (LP.Direction, x) -> Double -> [LP.Inequality x]
+constraintsFromSolution tol (dir,obj) opt =
+   case dir of
+      LP.Minimize -> [obj <=. opt + tol]
+      LP.Maximize -> [obj >=. opt - tol]
+
+successiveObjectives ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
+   Array sh a -> Double ->
+   NonEmpty.T [] (LP.Direction, [Term ix]) ->
+   ((LP.Direction, LP.Objective sh),
+    [(Double -> 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
+      (\(dir,obj) y1 ->
+         (constraintsFromSolution tol (dir,obj),
+          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)
+
+
+
+checkBound :: Double -> LP.Bound -> Double -> QC.Property
+checkBound tol bound x =
+   QC.counterexample (show (x, bound)) $
+   case bound of
+      LP.LessEqual up -> x<=up+tol
+      LP.GreaterEqual lo -> x>=lo-tol
+      LP.Between lo up -> lo-tol<=x && x<=up+tol
+      LP.Equal y -> approxReal tol x y
+      LP.Free -> True
+
+checkBounds ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
+   Double -> LP.Bounds ix -> Array sh Double -> QC.Property
+checkBounds tol bounds sol =
+   QC.conjoin $ map (\(ix,bnd) -> checkBound tol bnd (sol!ix)) $
+   BoxedArray.toAssociations $
+   BoxedArray.fromAssociations (LP.GreaterEqual 0) (Array.shape sol) $
+   map (\(LP.Inequality ix bnd) -> (ix,bnd)) bounds
+
+checkContraint ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
+   Double -> LP.Inequality [LP.Term Double ix] -> Array sh Double -> QC.Property
+checkContraint tol (LP.Inequality terms bnd) sol =
+   checkBound tol bnd $
+   scalarProductTerms (map (\(LP.Term c ix) -> (c,ix)) terms) sol
+
+checkFeasibility ::
+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
+   Double -> LP.Bounds ix -> Constraints ix -> Array sh Double -> QC.Property
+checkFeasibility tol bounds constrs sol =
+   checkBounds tol bounds sol
+   .&&.
+   QC.conjoin (map (flip (checkContraint tol) sol) constrs)
+
+
+affineCombination ::
+   (Shape.C sh, Eq sh, Storable a, Num a) =>
+   a -> Array sh a -> Array sh a -> Array sh a
+affineCombination c x y =
+   Array.zipWith (+) (Array.map ((1-c)*) x) (Array.map (c*) y)
+
+scalarProduct ::
+   (Shape.C sh, Eq sh, Storable a, Num a) =>
+   Array sh a -> Array sh a -> a
+scalarProduct x y = Array.sum $ Array.zipWith (*) x y