diff --git a/cbits/support.cpp b/cbits/support.cpp
--- a/cbits/support.cpp
+++ b/cbits/support.cpp
@@ -1,5 +1,5 @@
 #define CLP_EXTERN_C
-#include "coin/ClpSimplex.hpp"
+#include "ClpSimplex.hpp"
 #include "support.h"
 
 
diff --git a/coinor-clp.cabal b/coinor-clp.cabal
--- a/coinor-clp.cabal
+++ b/coinor-clp.cabal
@@ -1,6 +1,6 @@
 Cabal-Version:    2.2
 Name:             coinor-clp
-Version:          0.0.0.1
+Version:          0.0.0.2
 License:          BSD-3-Clause
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -35,6 +35,12 @@
 Extra-Source-Files:
   Makefile
   include/support.h
+  nix/clp-extern-c.patch
+  nix/coin-utils.nix
+  nix/osi.nix
+  nix/clp.nix
+  shell.nix
+  tool/ltrace.config
 
 Flag debug
   Description: Enable debug output
@@ -42,7 +48,7 @@
   Manual:      True
 
 Source-Repository this
-  Tag:         0.0.0.1
+  Tag:         0.0.0.2
   Type:        darcs
   Location:    https://hub.darcs.net/thielema/coinor-clp/
 
@@ -76,7 +82,7 @@
     Numeric.COINOR.CLP.Private
     Numeric.COINOR.CLP.Debug
 
-  PkgConfig-Depends: clp
+  PkgConfig-Depends: clp, coinutils
   Include-Dirs: include
   Extra-Libraries: stdc++
   Cxx-Sources:
@@ -86,7 +92,7 @@
   Type: exitcode-stdio-1.0
   Build-Depends:
     coinor-clp,
-    linear-programming,
+    linear-programming >= 0.0.1,
     comfort-array >=0.5.2,
     transformers,
     non-empty,
diff --git a/include/support.h b/include/support.h
--- a/include/support.h
+++ b/include/support.h
@@ -1,8 +1,8 @@
 #ifndef COINOR_HS_SUPPORT_H
 #define COINOR_HS_SUPPORT_H
 
-#include "coin/Clp_C_Interface.h"
-#include "coin/ClpPlusMinusOneMatrix.hpp"
+#include "Clp_C_Interface.h"
+#include "ClpPlusMinusOneMatrix.hpp"
 
 
 #ifdef __cplusplus
diff --git a/nix/clp-extern-c.patch b/nix/clp-extern-c.patch
new file mode 100644
--- /dev/null
+++ b/nix/clp-extern-c.patch
@@ -0,0 +1,23 @@
+diff --git a/src/Clp_C_Interface.h b/src/Clp_C_Interface.h
+index 567a8e28..005401db 100644
+--- a/src/Clp_C_Interface.h
++++ b/src/Clp_C_Interface.h
+@@ -18,8 +18,18 @@
+ #define CLP_LINKAGE_CB
+ #endif
+ 
++/** User does not need to see structure of model but C++ code does */
++#if defined(CLP_EXTERN_C)
++/* Real typedef for structure */
++class CMessageHandler;
++typedef struct {
++  ClpSimplex *model_;
++  CMessageHandler *handler_;
++} Clp_Simplex;
++#else
+ struct Clp_Simplex_s;
+ typedef struct Clp_Simplex_s Clp_Simplex;
++#endif
+ 
+ struct Clp_Solve_s;
+ typedef struct Clp_Solve_s Clp_Solve;
diff --git a/nix/clp.nix b/nix/clp.nix
new file mode 100644
--- /dev/null
+++ b/nix/clp.nix
@@ -0,0 +1,28 @@
+{ lib, stdenv, fetchFromGitHub, pkg-config, coin-utils, zlib, osi }:
+
+stdenv.mkDerivation rec {
+  version = "1.17.9-latest";
+  pname = "clp";
+  src = fetchFromGitHub {
+    owner = "coin-or";
+    repo = "Clp";
+    rev = "c9062fc02ccac08248cbf7b3ffae02994d7ea8c2";
+    hash = "sha256-vQ4X+8O5FyXjlHvxlfS6E1vA29ssvGes3pLkg5Zbdfo=";
+  };
+
+  patches = [ ./clp-extern-c.patch ];
+
+  nativeBuildInputs = [ pkg-config ];
+
+  propagatedBuildInputs = [ zlib coin-utils osi ];
+
+  doCheck = true;
+
+  meta = with lib; {
+    license = licenses.epl20;
+    homepage = "https://github.com/coin-or/Clp";
+    description = "An open-source linear programming solver written in C++";
+    platforms = platforms.darwin ++ platforms.linux;
+    maintainers = [ maintainers.vbgl ];
+  };
+}
diff --git a/nix/coin-utils.nix b/nix/coin-utils.nix
new file mode 100644
--- /dev/null
+++ b/nix/coin-utils.nix
@@ -0,0 +1,24 @@
+{ lib, stdenv, fetchFromGitHub, pkg-config, lapack }:
+
+stdenv.mkDerivation rec {
+  version = "2.11.10-latest";
+  pname = "coinutils";
+
+  src = fetchFromGitHub {
+    owner = "coin-or";
+    repo = "CoinUtils";
+    rev = "293e6e981774ed047e8f00f7aff9252262f83a02";
+    hash = "sha256-H3GBO5voOJO/0gMC0QY3me1qQFsCcOg5SuuoI3YV0IE=";
+  };
+
+  nativeBuildInputs = [ pkg-config lapack ];
+
+  doCheck = true;
+
+  meta = with lib; {
+    license = licenses.epl20;
+    homepage = "https://github.com/coin-or/CoinUtils";
+    description = "Collection of classes and helper functions that are generally useful to multiple COIN-OR projects";
+    maintainers = with maintainers; [ tmarkus ];
+  };
+}
diff --git a/nix/osi.nix b/nix/osi.nix
new file mode 100644
--- /dev/null
+++ b/nix/osi.nix
@@ -0,0 +1,44 @@
+{ stdenv, lib, fetchFromGitHub, gfortran, pkg-config
+, blas, zlib, bzip2, coin-utils
+, withGurobi ? false, gurobi
+, withCplex ? false, cplex }:
+
+stdenv.mkDerivation rec {
+  pname = "osi";
+  version = "0.108.9-latest";
+
+  src = fetchFromGitHub {
+    owner = "coin-or";
+    repo = "Osi";
+    rev = "477c9cd4f947936655feb3a5a3e38011052c1c25";
+    hash = "sha256-1yuDQRQth97ygNjiNtz9FwQO9584SO0iksFR0G9yNks=";
+  };
+
+  buildInputs =
+    [ blas zlib bzip2 coin-utils ]
+    ++ lib.optional withGurobi gurobi
+    ++ lib.optional withCplex cplex;
+  nativeBuildInputs = [ gfortran pkg-config ];
+  configureFlags =
+    lib.optionals withGurobi [ "--with-gurobi-incdir=${gurobi}/include" "--with-gurobi-lib=-lgurobi${gurobi.libSuffix}" ]
+    ++ lib.optionals withCplex [ "--with-cplex-incdir=${cplex}/cplex/include/ilcplex" "--with-cplex-lib=-lcplex${cplex.libSuffix}" ];
+
+  NIX_LDFLAGS =
+    lib.optionalString withCplex "-L${cplex}/cplex/bin/${cplex.libArch}";
+
+  # Compile errors
+  env.NIX_CFLAGS_COMPILE = "-Wno-cast-qual";
+  hardeningDisable = [ "format" ];
+
+  enableParallelBuilding = true;
+
+  passthru = { inherit withGurobi withCplex; };
+
+  meta = with lib; {
+    description = "An abstract base class to a generic linear programming (LP) solver";
+    homepage = "https://github.com/coin-or/Osi";
+    license = licenses.epl20;
+    platforms = platforms.unix;
+    maintainers = with maintainers; [ abbradar ];
+  };
+}
diff --git a/shell.nix b/shell.nix
new file mode 100644
--- /dev/null
+++ b/shell.nix
@@ -0,0 +1,34 @@
+let
+  pkgs = import <nixpkgs> {};
+  latest-coin-utils = pkgs.callPackage ./nix/coin-utils.nix {};
+  latest-osi =
+    (pkgs.callPackage ./nix/osi.nix {}).override {
+      coin-utils = latest-coin-utils;
+    };
+  latest-clp =
+    (pkgs.callPackage ./nix/clp.nix {}).override {
+      coin-utils = latest-coin-utils;
+      osi = latest-osi;
+    };
+in
+pkgs.mkShell {
+  # "override" required dependencies
+  buildInputs = (with pkgs; [
+    latest-clp
+    latest-coin-utils
+    pkg-config
+    gmp
+    bzip2
+    lapack
+    blas
+  ]);
+  shellHook = ''
+    # Why are these paths not set by Nix?
+    export PKG_CONFIG_PATH=${latest-coin-utils}/lib/pkgconfig:$PKG_CONFIG_PATH
+    export LD_LIBRARY_PATH=${latest-coin-utils}/lib:$LD_LIBRARY_PATH
+    export PKG_CONFIG_PATH=${latest-osi}/lib/pkgconfig:$PKG_CONFIG_PATH
+    export LD_LIBRARY_PATH=${latest-osi}/lib:$LD_LIBRARY_PATH
+    export PKG_CONFIG_PATH=${latest-clp}/lib/pkgconfig:$PKG_CONFIG_PATH
+    export LD_LIBRARY_PATH=${latest-clp}/lib:$LD_LIBRARY_PATH
+  '';
+}
diff --git a/src/Numeric/COINOR/CLP.hs b/src/Numeric/COINOR/CLP.hs
--- a/src/Numeric/COINOR/CLP.hs
+++ b/src/Numeric/COINOR/CLP.hs
@@ -21,7 +21,8 @@
          (Method(runMethod), Result, FailureType(..),
           runContT, withBuffer, false,
           storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,
-          storeConstraints, setOptimizationDirection, examineStatus)
+          storeConstraints, prepareConstraints,
+          setOptimizationDirection, examineStatus)
 
 import qualified Numeric.LinearProgramming.Common as LP
 import Numeric.LinearProgramming.Common
@@ -45,11 +46,13 @@
 {- $setup
 >>> import qualified Numeric.COINOR.CLP as LP
 >>> import qualified Numeric.LinearProgramming.Test as TestLP
->>> import Numeric.COINOR.CLP (PlusMinusOne(..), (.*), (==.), (<=.), (>=<.))
+>>> import Numeric.COINOR.CLP
+>>>    (PlusMinusOne(..), (.*), (==.), (<=.), (>=.), (>=<.))
 >>>
 >>> import qualified Data.Array.Comfort.Storable as Array
 >>> import qualified Data.Array.Comfort.Shape as Shape
 >>>
+>>> import Data.Either.HT (mapRight)
 >>> import Data.Tuple.HT (mapSnd)
 >>>
 >>> import qualified Test.QuickCheck as QC
@@ -103,14 +106,15 @@
 
 instance Coefficient Double where
    loadProblem shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr = do
-      (coefficientsShape, coefficientsPtr, indexPtr, startPtr)
-         <- storeConstraints shape constrs
+      let (coefficients, indices, rowStarts) = prepareConstraints shape constrs
+      (coefficientsPtr, indexPtr, startPtr)
+         <- storeConstraints (coefficients, indices, rowStarts)
       let createMatrix =
             FFI.newCoinPackedMatrix
                false
                (fromIntegral $ Shape.size shape)
                (fromIntegral $ length constrs)
-               (fromIntegral $ Shape.size coefficientsShape)
+               (fromIntegral $ Shape.size $ Array.shape coefficients)
                coefficientsPtr
                indexPtr
                startPtr
@@ -171,18 +175,22 @@
          mapSnd Array.toTuple <$>
          LP.simplex LP.dual []
             [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20]
-            (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
+            (LP.Maximize, Array.fromTuple (4,-3,2)
+               :: Array.Array TripletShape Double)
 :}
 Right (28.0,(5.0,0.0,4.0))
+
 >>> :{
    case Shape.indexTupleFromShape tripletShape of
       (x,y,z) ->
          mapSnd Array.toTuple <$>
          LP.simplex LP.primal [y >=<. (-12,12)]
             [[1.*x, (-1).*y] <=. 10, [(-1).*y, (1::Double).*z] <=. 20]
-            (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
+            (LP.Maximize, Array.fromTuple (4,-3,2)
+               :: Array.Array TripletShape Double)
 :}
 Right (116.0,(22.0,12.0,32.0))
+
 >>> :{
    case Shape.indexTupleFromShape tripletShape of
       (x,y,z) ->
@@ -190,10 +198,33 @@
          LP.simplex LP.primal [y >=<. (-12,12)]
             [[PlusOne .* x, MinusOne .* y] <=. 10,
              [MinusOne .* y, PlusOne .* z] <=. 20]
-            (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
+            (LP.Maximize, Array.fromTuple (4,-3,2)
+               :: Array.Array TripletShape Double)
 :}
 Right (116.0,(22.0,12.0,32.0))
 
+>>> :{
+   case Shape.indexTupleFromShape tripletShape of
+      (x,y,z) ->
+         mapSnd Array.toTuple <$>
+         LP.simplex LP.primal [y >=<. (-12,12)]
+            [[1.*x, 1.*y] <=. 10, [1.*y, (-1::Double).*z] >=. 20]
+            (LP.Maximize, Array.fromTuple (4,3,2)
+               :: Array.Array TripletShape Double)
+:}
+Left PrimalInfeasible
+
+>>> :{
+   case Shape.indexTupleFromShape tripletShape of
+      (x,y,z) ->
+         mapSnd Array.toTuple <$>
+         LP.simplex LP.primal [y >=<. (-12,12)]
+            [[1.*x, 1.*y] <=. 10, [1.*y, (1::Double).*z] >=. 20]
+            (LP.Maximize, Array.fromTuple (4,3,2)
+               :: Array.Array TripletShape Double)
+:}
+Left DualInfeasible
+
 prop> :{
    forAllMethod $ \method (QC.Positive posWeight) (QC.Positive negWeight) target ->
    case Shape.indexTupleFromShape pairShape of
@@ -270,6 +301,59 @@
             in case dir of
                   LP.Minimize -> opt-0.01 <= val
                   LP.Maximize -> opt+0.01 >= val
+:}
+prop> :{
+   forAllMethod $ \method ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllBoundedProblem origin $ \bounds constrs ->
+   QC.forAll (TestLP.genObjective origin) $ \dirObjA ->
+   QC.forAll (TestLP.genObjective origin) $ \dirObjB ->
+   let solA = LP.simplex method bounds constrs dirObjA in
+   let solB = LP.simplex method bounds constrs dirObjB in
+   QC.counterexample (show (mapRight fst solA, mapRight fst solB)) $
+   case (solA, solB) of
+      (Right _, Right _) -> True
+      (Left _, Left _) -> True
+      _ -> False
+:}
+prop> :{
+   forAllMethod $ \method ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllProblem origin $ \bounds constrs ->
+   QC.forAll (TestLP.genObjective origin) $ \(_dir,obj) ->
+   case (LP.simplex method bounds constrs (LP.Minimize,obj),
+         LP.simplex method bounds constrs (LP.Maximize,obj)) of
+      (Right (optMin,_), Right (optMax,_)) ->
+         QC.counterexample (show (optMin, optMax)) $ optMin <= optMax + 0.01
+      _ -> QC.property False
+:}
+prop> :{
+   forAllMethod $ \method ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllProblem origin $ \bounds allConstrs ->
+   QC.forAll (QC.sublistOf allConstrs) $ \someConstrs ->
+   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
+   case (LP.simplex method bounds allConstrs (dir,obj),
+         LP.simplex method bounds someConstrs (dir,obj)) of
+      (Right (optAll,_), Right (optSome,_)) ->
+         QC.counterexample (show (optAll, optSome)) $
+         case dir of
+            LP.Minimize -> optAll >= optSome-0.01
+            LP.Maximize -> optAll <= optSome+0.01
+      _ -> QC.property False
+:}
+prop> :{
+   forAllMethod $ \methodA ->
+   forAllMethod $ \methodB ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllProblem origin $ \bounds constrs ->
+   QC.forAll (TestLP.genObjective origin) $ \dirObj ->
+   case (LP.simplex methodA bounds constrs dirObj,
+         LP.simplex methodB bounds constrs dirObj) of
+      (Right (optA,_), Right (optB,_)) ->
+         QC.counterexample (show (optA, optB)) $
+         approxReal 0.01 optA optB
+      _ -> QC.property False
 :}
 -}
 simplex ::
diff --git a/src/Numeric/COINOR/CLP/FFI.hsc b/src/Numeric/COINOR/CLP/FFI.hsc
--- a/src/Numeric/COINOR/CLP/FFI.hsc
+++ b/src/Numeric/COINOR/CLP/FFI.hsc
@@ -2,6 +2,7 @@
 module Numeric.COINOR.CLP.FFI where
 
 import qualified Foreign.C.Types as C
+import Foreign.C.String (CString)
 import Foreign.Ptr (Ptr)
 
 import Data.Int (Int32)
@@ -12,7 +13,7 @@
 type CBool = C.CBool
 
 
-#include "coin/Clp_C_Interface.h"
+#include "Clp_C_Interface.h"
 
 type BigIndex = #{type CoinBigIndex}
 
@@ -157,3 +158,12 @@
 
 foreign import ccall "Clp_dumpMatrix"
    dumpMatrix :: Ptr Simplex -> IO ()
+
+
+newtype CoinorBool = CoinorBool CInt
+
+foreign import ccall "Clp_readMps"
+   readMps :: Ptr Simplex -> CString -> CoinorBool -> CoinorBool -> IO CInt
+
+foreign import ccall "Clp_writeMps"
+   writeMps :: Ptr Simplex -> CString -> CInt -> CInt -> CDouble -> IO CInt
diff --git a/src/Numeric/COINOR/CLP/Monad.hs b/src/Numeric/COINOR/CLP/Monad.hs
--- a/src/Numeric/COINOR/CLP/Monad.hs
+++ b/src/Numeric/COINOR/CLP/Monad.hs
@@ -20,10 +20,11 @@
          (Method(runMethod), Result,
           runContT, withBuffer,
           storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,
-          storeConstraints, setOptimizationDirection, examineStatus)
+          storeConstraints, prepareConstraints,
+          setOptimizationDirection, examineStatus)
 
 import Numeric.LinearProgramming.Common
-         (Bounds, Constraints, Direction(..), Objective, Term(Term))
+         (Bounds, Constraints, Direction(..), Objective)
 
 import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
@@ -173,7 +174,7 @@
 
 {- |
 Add new constraints to an existing problem
-and simplex with a new direction and objective.
+and run with a new direction and objective.
 
 >>> :{
    case Shape.indexTupleFromShape tripletShape of
@@ -231,15 +232,8 @@
       error "COINOR.CLP.Monad.solve: objective shape mismatch"
 
    liftIO $ runContT $ do
-      {- ToDo: Remove this in the future.
-      This is a work-around for a bug in Clp_addRows
-      https://hydra.nixos.org/build/239790474/nixlog/2
-      https://list.coin-or.org/pipermail/clp/2023-November/001805.html
-      -}
-      let constrsFiltered =
-            map (fmap (filter (\(Term c _x) -> c/=0)))constrs
-      (_, coefficientsPtr, indexPtr, startPtr) <-
-         storeConstraints shape constrsFiltered
+      (coefficientsPtr, indexPtr, startPtr) <-
+         storeConstraints $ prepareConstraints shape constrs
       (rowlbPtr,rowubPtr) <- storeBounds $ prepareRowBoundsArrays constrs
       objPtr <- withBuffer $ Array.map realToFrac obj
       liftIO $ do
diff --git a/src/Numeric/COINOR/CLP/Private.hs b/src/Numeric/COINOR/CLP/Private.hs
--- a/src/Numeric/COINOR/CLP/Private.hs
+++ b/src/Numeric/COINOR/CLP/Private.hs
@@ -80,29 +80,43 @@
 storeBounds = AppHT.mapPair (withBuffer, withBuffer)
 
 
-storeConstraints ::
+prepareConstraints ::
    (Shape.Indexed sh, Shape.Index sh ~ ix) =>
    sh -> Constraints Double ix ->
-   MC.ContT r IO (ShapeInt, Ptr CDouble, Ptr CInt, Ptr FFI.BigIndex)
-storeConstraints shape constrs = do
-   let rowStarts =
+   (Array ShapeInt CDouble, Array ShapeInt CInt, Array ShapeInt FFI.BigIndex)
+prepareConstraints shape constrs =
+   let {-
+       It seems that LP solvers generally do not expect zero coefficients,
+       although that is not document.
+       https://hydra.nixos.org/build/239790474/nixlog/2
+       https://list.coin-or.org/pipermail/clp/2023-November/001805.html
+       -}
+       constrsNonZero = map (fmap (filter (\(Term c _x) -> c/=0))) constrs
+       rowStarts =
          Array.vectorFromList $ scanl (+) 0 $
          map (\(Inequality terms _bnd) -> fromIntegral $ length terms)
-            constrs
-   startPtr <- withBuffer rowStarts
-   let shapeOffset = Shape.offset shape
-   let coefficients =
-         concatMap (\(Inequality terms _bnd) -> terms) constrs
-   let coefficientArr =
+            constrsNonZero
+       shapeOffset = Shape.offset shape
+       coefficients =
+         concatMap (\(Inequality terms _bnd) -> terms) constrsNonZero
+       indexArr =
          Array.vectorFromList $
          map (\(Term _ ix) -> fromIntegral $ shapeOffset ix) coefficients
-   indexPtr <- withBuffer coefficientArr
-   coefficientsPtr <-
-      withBuffer $ Array.vectorFromList $
-      map (\(Term c _) -> realToFrac c) coefficients
-   return (Array.shape coefficientArr, coefficientsPtr, indexPtr, startPtr)
+       coefficientArr =
+         Array.vectorFromList $
+         map (\(Term c _) -> realToFrac c) coefficients
+   in (coefficientArr, indexArr, rowStarts)
 
+storeConstraints ::
+   (Array ShapeInt CDouble, Array ShapeInt CInt, Array ShapeInt FFI.BigIndex) ->
+   MC.ContT r IO (Ptr CDouble, Ptr CInt, Ptr FFI.BigIndex)
+storeConstraints (coefficients, indices, rowStarts) =
+   AppHT.lift3 (,,)
+      (withBuffer coefficients)
+      (withBuffer indices)
+      (withBuffer rowStarts)
 
+
 setOptimizationDirection :: Ptr FFI.Simplex -> Direction -> IO ()
 setOptimizationDirection lp dir =
    FFI.setOptimizationDirection lp $
@@ -133,9 +147,7 @@
 
 type Result sh = Either FailureType (Double, Array sh Double)
 
-examineStatus ::
-   (Shape.C sh) =>
-   sh -> Ptr FFI.Simplex -> IO (Either FailureType (Double, Array sh Double))
+examineStatus :: (Shape.C sh) => sh -> Ptr FFI.Simplex -> IO (Result sh)
 examineStatus shape lp = do
    status <- FFI.status lp
    case status of
diff --git a/test/Test/Numeric/COINOR/CLP.hs b/test/Test/Numeric/COINOR/CLP.hs
--- a/test/Test/Numeric/COINOR/CLP.hs
+++ b/test/Test/Numeric/COINOR/CLP.hs
@@ -1,19 +1,21 @@
 -- Do not edit! Automatically created with doctest-extract from src/Numeric/COINOR/CLP.hs
-{-# LINE 45 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 46 "src/Numeric/COINOR/CLP.hs" #-}
 
 module Test.Numeric.COINOR.CLP where
 
 import Test.DocTest.Base
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 46 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 47 "src/Numeric/COINOR/CLP.hs" #-}
 import     qualified Numeric.COINOR.CLP as LP
 import     qualified Numeric.LinearProgramming.Test as TestLP
-import     Numeric.COINOR.CLP (PlusMinusOne(..), (.*), (==.), (<=.), (>=<.))
+import     Numeric.COINOR.CLP
+       (PlusMinusOne(..), (.*), (==.), (<=.), (>=.), (>=<.))
 
 import     qualified Data.Array.Comfort.Storable as Array
 import     qualified Data.Array.Comfort.Shape as Shape
 
+import     Data.Either.HT (mapRight)
 import     Data.Tuple.HT (mapSnd)
 
 import     qualified Test.QuickCheck as QC
@@ -50,36 +52,38 @@
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Numeric.COINOR.CLP:168: "
-{-# LINE 168 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:172: "
+{-# LINE 172 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.example(
-{-# LINE 168 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 172 "src/Numeric/COINOR/CLP.hs" #-}
       
    case Shape.indexTupleFromShape tripletShape of
       (x,y,z) ->
          mapSnd Array.toTuple <$>
          LP.simplex LP.dual []
             [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20]
-            (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
+            (LP.Maximize, Array.fromTuple (4,-3,2)
+               :: Array.Array TripletShape Double)
   )
   [ExpectedLine [LineChunk "Right (28.0,(5.0,0.0,4.0))"]]
- DocTest.printPrefix "Numeric.COINOR.CLP:177: "
-{-# LINE 177 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:183: "
+{-# LINE 183 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.example(
-{-# LINE 177 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 183 "src/Numeric/COINOR/CLP.hs" #-}
       
    case Shape.indexTupleFromShape tripletShape of
       (x,y,z) ->
          mapSnd Array.toTuple <$>
          LP.simplex LP.primal [y >=<. (-12,12)]
             [[1.*x, (-1).*y] <=. 10, [(-1).*y, (1::Double).*z] <=. 20]
-            (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
+            (LP.Maximize, Array.fromTuple (4,-3,2)
+               :: Array.Array TripletShape Double)
   )
   [ExpectedLine [LineChunk "Right (116.0,(22.0,12.0,32.0))"]]
- DocTest.printPrefix "Numeric.COINOR.CLP:186: "
-{-# LINE 186 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:194: "
+{-# LINE 194 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.example(
-{-# LINE 186 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 194 "src/Numeric/COINOR/CLP.hs" #-}
       
    case Shape.indexTupleFromShape tripletShape of
       (x,y,z) ->
@@ -87,13 +91,42 @@
          LP.simplex LP.primal [y >=<. (-12,12)]
             [[PlusOne .* x, MinusOne .* y] <=. 10,
              [MinusOne .* y, PlusOne .* z] <=. 20]
-            (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)
+            (LP.Maximize, Array.fromTuple (4,-3,2)
+               :: Array.Array TripletShape Double)
   )
   [ExpectedLine [LineChunk "Right (116.0,(22.0,12.0,32.0))"]]
- DocTest.printPrefix "Numeric.COINOR.CLP:197: "
-{-# LINE 197 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:206: "
+{-# LINE 206 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.example(
+{-# LINE 206 "src/Numeric/COINOR/CLP.hs" #-}
+      
+   case Shape.indexTupleFromShape tripletShape of
+      (x,y,z) ->
+         mapSnd Array.toTuple <$>
+         LP.simplex LP.primal [y >=<. (-12,12)]
+            [[1.*x, 1.*y] <=. 10, [1.*y, (-1::Double).*z] >=. 20]
+            (LP.Maximize, Array.fromTuple (4,3,2)
+               :: Array.Array TripletShape Double)
+  )
+  [ExpectedLine [LineChunk "Left PrimalInfeasible"]]
+ DocTest.printPrefix "Numeric.COINOR.CLP:217: "
+{-# LINE 217 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.example(
+{-# LINE 217 "src/Numeric/COINOR/CLP.hs" #-}
+      
+   case Shape.indexTupleFromShape tripletShape of
+      (x,y,z) ->
+         mapSnd Array.toTuple <$>
+         LP.simplex LP.primal [y >=<. (-12,12)]
+            [[1.*x, 1.*y] <=. 10, [1.*y, (1::Double).*z] >=. 20]
+            (LP.Maximize, Array.fromTuple (4,3,2)
+               :: Array.Array TripletShape Double)
+  )
+  [ExpectedLine [LineChunk "Left DualInfeasible"]]
+ DocTest.printPrefix "Numeric.COINOR.CLP:228: "
+{-# LINE 228 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.property(
-{-# LINE 197 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 228 "src/Numeric/COINOR/CLP.hs" #-}
         
    forAllMethod $ \method (QC.Positive posWeight) (QC.Positive negWeight) target ->
    case Shape.indexTupleFromShape pairShape of
@@ -109,10 +142,10 @@
                .&&.
                (posResult === 0 .||. negResult === 0)
   )
- DocTest.printPrefix "Numeric.COINOR.CLP:212: "
-{-# LINE 212 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:243: "
+{-# LINE 243 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.property(
-{-# LINE 212 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 243 "src/Numeric/COINOR/CLP.hs" #-}
         
    forAllMethod $ \method target ->
    case Shape.indexTupleFromShape pairShape of
@@ -129,10 +162,10 @@
                .&&.
                (posResult === 0 .||. negResult === 0)
   )
- DocTest.printPrefix "Numeric.COINOR.CLP:229: "
-{-# LINE 229 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:260: "
+{-# LINE 260 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.property(
-{-# LINE 229 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 260 "src/Numeric/COINOR/CLP.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
@@ -142,10 +175,10 @@
       Left _ -> False
       Right _ -> True
   )
- DocTest.printPrefix "Numeric.COINOR.CLP:238: "
-{-# LINE 238 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:269: "
+{-# LINE 269 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.property(
-{-# LINE 238 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 269 "src/Numeric/COINOR/CLP.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
@@ -155,10 +188,10 @@
       Left _ -> QC.property False
       Right (_,sol) -> TestLP.checkFeasibility 0.1 bounds constrs sol
   )
- DocTest.printPrefix "Numeric.COINOR.CLP:247: "
-{-# LINE 247 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:278: "
+{-# LINE 278 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.property(
-{-# LINE 247 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 278 "src/Numeric/COINOR/CLP.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
@@ -171,10 +204,10 @@
          TestLP.checkFeasibility 0.1 bounds constrs $
          TestLP.affineCombination lambda sol (Array.map fromIntegral origin)
   )
- DocTest.printPrefix "Numeric.COINOR.CLP:259: "
-{-# LINE 259 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP:290: "
+{-# LINE 290 "src/Numeric/COINOR/CLP.hs" #-}
  DocTest.property(
-{-# LINE 259 "src/Numeric/COINOR/CLP.hs" #-}
+{-# LINE 290 "src/Numeric/COINOR/CLP.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
@@ -189,4 +222,73 @@
             in case dir of
                   LP.Minimize -> opt-0.01 <= val
                   LP.Maximize -> opt+0.01 >= val
+  )
+ DocTest.printPrefix "Numeric.COINOR.CLP:305: "
+{-# LINE 305 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.property(
+{-# LINE 305 "src/Numeric/COINOR/CLP.hs" #-}
+        
+   forAllMethod $ \method ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllBoundedProblem origin $ \bounds constrs ->
+   QC.forAll (TestLP.genObjective origin) $ \dirObjA ->
+   QC.forAll (TestLP.genObjective origin) $ \dirObjB ->
+   let solA = LP.simplex method bounds constrs dirObjA in
+   let solB = LP.simplex method bounds constrs dirObjB in
+   QC.counterexample (show (mapRight fst solA, mapRight fst solB)) $
+   case (solA, solB) of
+      (Right _, Right _) -> True
+      (Left _, Left _) -> True
+      _ -> False
+  )
+ DocTest.printPrefix "Numeric.COINOR.CLP:319: "
+{-# LINE 319 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.property(
+{-# LINE 319 "src/Numeric/COINOR/CLP.hs" #-}
+        
+   forAllMethod $ \method ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllProblem origin $ \bounds constrs ->
+   QC.forAll (TestLP.genObjective origin) $ \(_dir,obj) ->
+   case (LP.simplex method bounds constrs (LP.Minimize,obj),
+         LP.simplex method bounds constrs (LP.Maximize,obj)) of
+      (Right (optMin,_), Right (optMax,_)) ->
+         QC.counterexample (show (optMin, optMax)) $ optMin <= optMax + 0.01
+      _ -> QC.property False
+  )
+ DocTest.printPrefix "Numeric.COINOR.CLP:330: "
+{-# LINE 330 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.property(
+{-# LINE 330 "src/Numeric/COINOR/CLP.hs" #-}
+        
+   forAllMethod $ \method ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllProblem origin $ \bounds allConstrs ->
+   QC.forAll (QC.sublistOf allConstrs) $ \someConstrs ->
+   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
+   case (LP.simplex method bounds allConstrs (dir,obj),
+         LP.simplex method bounds someConstrs (dir,obj)) of
+      (Right (optAll,_), Right (optSome,_)) ->
+         QC.counterexample (show (optAll, optSome)) $
+         case dir of
+            LP.Minimize -> optAll >= optSome-0.01
+            LP.Maximize -> optAll <= optSome+0.01
+      _ -> QC.property False
+  )
+ DocTest.printPrefix "Numeric.COINOR.CLP:345: "
+{-# LINE 345 "src/Numeric/COINOR/CLP.hs" #-}
+ DocTest.property(
+{-# LINE 345 "src/Numeric/COINOR/CLP.hs" #-}
+        
+   forAllMethod $ \methodA ->
+   forAllMethod $ \methodB ->
+   TestLP.forAllOrigin $ \origin ->
+   TestLP.forAllProblem origin $ \bounds constrs ->
+   QC.forAll (TestLP.genObjective origin) $ \dirObj ->
+   case (LP.simplex methodA bounds constrs dirObj,
+         LP.simplex methodB bounds constrs dirObj) of
+      (Right (optA,_), Right (optB,_)) ->
+         QC.counterexample (show (optA, optB)) $
+         approxReal 0.01 optA optB
+      _ -> QC.property False
   )
diff --git a/test/Test/Numeric/COINOR/CLP/Monad.hs b/test/Test/Numeric/COINOR/CLP/Monad.hs
--- a/test/Test/Numeric/COINOR/CLP/Monad.hs
+++ b/test/Test/Numeric/COINOR/CLP/Monad.hs
@@ -1,5 +1,5 @@
 -- Do not edit! Automatically created with doctest-extract from src/Numeric/COINOR/CLP/Monad.hs
-{-# LINE 42 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+{-# LINE 43 "src/Numeric/COINOR/CLP/Monad.hs" #-}
 
 {-# OPTIONS_GHC -XTypeFamilies #-}
 {-# OPTIONS_GHC -XTypeOperators #-}
@@ -8,7 +8,7 @@
 import Test.DocTest.Base
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 45 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+{-# LINE 46 "src/Numeric/COINOR/CLP/Monad.hs" #-}
 import     qualified Numeric.COINOR.CLP.Monad as LP
 import     qualified Numeric.COINOR.CLP as CLP
 import     Test.Numeric.COINOR.CLP.Utility (traverse_Lag, traverseLag)
@@ -114,10 +114,10 @@
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Numeric.COINOR.CLP.Monad:178: "
-{-# LINE 178 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:179: "
+{-# LINE 179 "src/Numeric/COINOR/CLP/Monad.hs" #-}
  DocTest.example(
-{-# LINE 178 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+{-# LINE 179 "src/Numeric/COINOR/CLP/Monad.hs" #-}
       
    case Shape.indexTupleFromShape tripletShape of
       (x,y,z) ->
@@ -128,10 +128,10 @@
                (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))
   )
   [ExpectedLine [LineChunk "Right (28.0,(5.0,0.0,4.0))"]]
- DocTest.printPrefix "Numeric.COINOR.CLP.Monad:189: "
-{-# LINE 189 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:190: "
+{-# LINE 190 "src/Numeric/COINOR/CLP/Monad.hs" #-}
  DocTest.property(
-{-# LINE 189 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+{-# LINE 190 "src/Numeric/COINOR/CLP/Monad.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
@@ -143,10 +143,10 @@
       (Right (optA,_), Right (optB,_)) ->
          TestLP.approxReal 0.1 optA optB; _ -> False
   )
- DocTest.printPrefix "Numeric.COINOR.CLP.Monad:201: "
-{-# LINE 201 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:202: "
+{-# LINE 202 "src/Numeric/COINOR/CLP/Monad.hs" #-}
  DocTest.property(
-{-# LINE 201 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+{-# LINE 202 "src/Numeric/COINOR/CLP/Monad.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
@@ -157,10 +157,10 @@
          either (\msg -> QC.counterexample (show msg) False) (const $ QC.property True) $
          runSuccessive method (Array.shape origin) bounds (constrs,dirObj) objs
   )
- DocTest.printPrefix "Numeric.COINOR.CLP.Monad:212: "
-{-# LINE 212 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+ DocTest.printPrefix "Numeric.COINOR.CLP.Monad:213: "
+{-# LINE 213 "src/Numeric/COINOR/CLP/Monad.hs" #-}
  DocTest.property(
-{-# LINE 212 "src/Numeric/COINOR/CLP/Monad.hs" #-}
+{-# LINE 213 "src/Numeric/COINOR/CLP/Monad.hs" #-}
         
    forAllMethod $ \method ->
    TestLP.forAllOrigin $ \origin ->
diff --git a/tool/ltrace.config b/tool/ltrace.config
new file mode 100644
--- /dev/null
+++ b/tool/ltrace.config
@@ -0,0 +1,25 @@
+typedef Clp_Simplex = addr
+typedef CoinBigIndex = int
+
+Clp_Simplex Clp_newModel();
+void Clp_deleteModel(Clp_Simplex);
+
+void Clp_setOptimizationDirection(Clp_Simplex, double);
+
+void Clp_loadProblem(Clp_Simplex, int, int, array(int,arg3)*, int*, double*, array(double,arg2)*, array(double,arg2)*, array(double,arg2)*, array(double,arg3)*, array(double,arg3)*);
+
+
+void Clp_addColumns(Clp_Simplex, int, array(double,arg2)*, array(double,arg2)*, double *, CoinBigIndex *, int *, double *);
+void Clp_addRows(Clp_Simplex, int, array(double,arg2)*, array(double,arg2)*, CoinBigIndex *, int *, double *);
+
+void Clp_chgObjCoefficients(Clp_Simplex, double *);
+
+int Clp_dual(Clp_Simplex, int);
+int Clp_primal(Clp_Simplex, int);
+
+double Clp_objectiveValue(Clp_Simplex);
+double *Clp_getColSolution(Clp_Simplex);
+
+void Clp_setLogLevel(Clp_Simplex, int);
+int Clp_status(Clp_Simplex);
+void Clp_dumpMatrix (Clp_Simplex);