2 files changed, 27 insertions, 111 deletions

diff --git a/packages/glpk/lib/Numeric/LinearProgramming.hs b/packages/glpk/lib/Numeric/LinearProgramming.hs
index a888214..78f55a2 100644
--- a/packages/glpk/lib/Numeric/LinearProgramming.hs
+++ b/packages/glpk/lib/Numeric/LinearProgramming.hs
@@ -28,8 +28,8 @@ can be solved as follows:
 
 prob = Maximize [4, -3, 2]
 
-constr1 = Sparse [ [2#1, 1#2] :<: 10
-                 , [1#2, 5#3] :<: 20
+constr1 = Sparse [ [2\#1, 1\#2] :<: 10
+                 , [1\#2, 5\#3] :<: 20
                  ]
 
 \> simplex prob constr1 []
@@ -52,6 +52,8 @@ Unbounded@
 
 The given bound for a variable completely replaces the default,
 so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@.
+Multiple bounds for a variable are not allowed, instead of
+@[i :>: a, i:<: b]@ use @i :&: (a,b)@.
 
 -}
 
@@ -65,11 +67,11 @@ module Numeric.LinearProgramming(
     Solution(..)
 ) where
 
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra hiding (i)
 import Data.Packed.Development
 import Foreign(Ptr,unsafePerformIO)
 import Foreign.C.Types(CInt)
-import Data.List((\\),sortBy)
+import Data.List((\\),sortBy,nub)
 import Data.Function(on)
 
 --import Debug.Trace
@@ -82,7 +84,7 @@ import Data.Function(on)
 infixl 5 #
 (#) = (,)
 
-data Bound x =  x :<: Double 
+data Bound x =  x :<: Double
              |  x :>: Double
              |  x :&: (Double,Double)
              |  x :==: Double
@@ -107,11 +109,13 @@ type Bounds = [Bound Int]
 
 simplex :: Optimization -> Constraints -> Bounds -> Solution
 
-simplex opt (Dense  []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+simplex opt (Dense  []) bnds = simplex opt (Sparse []) bnds
 simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
 
 simplex opt (Dense constr) bnds = extract sg sol where
-    sol = simplexDense (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
+    sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
+    n = length objfun
+    m = length constr
     (sz, sg, objfun) = adapt opt
 
 simplex opt (Sparse constr) bnds = extract sg sol where
@@ -178,46 +182,37 @@ mkBound2 b = (obj b, mkBound1 b)
 
 mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double
 mkBounds n b1 b2 = fromLists (cb++vb) where
-    gv = map obj b2
+    gv' = map obj b2
+    gv | nub gv' == gv' = gv'
+       | otherwise = error $ "simplex: duplicate bounds for vars " ++ show (gv'\\nub gv')
     rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv
        | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n
     vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2
     cb = map mkBound1 b1
 
 mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double
-mkConstrD n f b1 | ok = fromLists (f : cs)
+mkConstrD n f b1 | ok = fromLists (ob ++ co)
                  | otherwise = error $ "simplex: dense constraints require "++show n
                                      ++" variables, given " ++ show ls
     where
-       cs | null b1   = error $ "simplex: dense: empty constraints"
-          | otherwise = map obj b1
+       cs = map obj b1
        ls = map length cs
        ok = all (==n) ls
+       den = fromLists cs
+       ob = map (([0,0]++).return) f
+       co = [[fromIntegral i, fromIntegral j,den@@>(i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]
 
 mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double
 mkConstrS n objfun b1 = fromLists (ob ++ co) where
     ob = map (([0,0]++).return) objfun
     co = concat $ zipWith f [1::Int ..] cs
-    cs    | null b1   = error $ "simplex: sparse: empty constraints"
-          | otherwise = map obj b1
+    cs = map obj b1
     f k = map (g k)
     g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c]
               | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n
 
 -----------------------------------------------------
 
-foreign import ccall "c_simplex_dense" c_simplex_dense
-    :: CInt -> CInt -> Ptr Double    -- coeffs
-    -> CInt -> CInt -> Ptr Double    -- bounds
-    -> CInt -> Ptr Double            -- result
-    -> IO CInt                       -- exit code
-
-simplexDense :: Matrix Double -> Matrix Double -> Vector Double
-simplexDense c b = unsafePerformIO $ do
-    s <- createVector (2+cols c)
-    app3 c_simplex_dense mat (cmat c) mat (cmat b) vec s "c_simplex_dense"
-    return s
-
 foreign import ccall "c_simplex_sparse" c_simplex_sparse
     :: CInt -> CInt                  -- rows and cols
     -> CInt -> CInt -> Ptr Double    -- coeffs
@@ -241,14 +236,6 @@ glpFX = 4
 
 {- Raw format of coeffs
 
-simplexDense
-
-((3+1) >< 3)
-  [ 10, 6, 4
-  ,  1, 1, 1
-  , 10, 4, 5
-  ,  2, 2, 6 :: Double]
-
 simplexSparse
 
 (12><3)
@@ -264,7 +251,7 @@ simplexSparse
  , 3.0, 1.0,  2.0
  , 3.0, 2.0,  2.0
  , 3.0, 3.0,  6.0 ]
-  
+
 bounds = (6><3)
   [ glpUP,0,100
   , glpUP,0,600
diff --git a/packages/glpk/lib/Numeric/LinearProgramming/glpk.c b/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
index 0f4ad79..bfbb435 100644
--- a/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
+++ b/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
@@ -1,84 +1,14 @@
-typedef struct { double r, i; } doublecomplex;
-
 #define DVEC(A) int A##n, double*A##p
-#define CVEC(A) int A##n, doublecomplex*A##p
 #define DMAT(A) int A##r, int A##c, double*A##p
-#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
 
 #define AT(M,r,co) (M##p[(r)*M##c+(co)])
 
-/*-----------------------------------------------------*/
-
 #include <stdlib.h>
 #include <stdio.h>
 #include <glpk.h>
 #include <math.h>
 
-int c_simplex_dense(DMAT(c), DMAT(b), DVEC(s)) {
-    glp_prob *lp;
-    lp = glp_create_prob();
-    glp_set_obj_dir(lp, GLP_MAX);
-
-    int m = cr-1;
-    int n = cc;
-    glp_add_rows(lp, m);
-    glp_add_cols(lp, n);
-    int * ia = malloc((1+m*n)*sizeof(int));
-    int * ja = malloc((1+m*n)*sizeof(int));
-    double * ar = malloc((1+m*n)*sizeof(double));
-    int i,j,k;
-    k=0;
-    for (i=1;i<=m;i++) {
-        for(j=1;j<=n;j++) {
-            k++;
-            ia[k] = i;
-            ja[k] = j;
-            ar[k] = AT(c,i,j-1);
-            //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
-        }
-    }
-    glp_load_matrix(lp, k, ia, ja, ar);
-    for (j=1;j<=n;j++) {
-        glp_set_obj_coef(lp, j, AT(c,0,j-1));
-    }
-    int t;
-    for (i=1;i<=m;i++) {
-    switch((int)rint(AT(b,i-1,0))) {
-        case 0: { t = GLP_FR; break; } 
-        case 1: { t = GLP_LO; break; }
-        case 2: { t = GLP_UP; break; }
-        case 3: { t = GLP_DB; break; }
-       default: { t = GLP_FX; break; }  
-    }
-    glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
-    }
-    for (j=1;j<=n;j++) {
-    switch((int)rint(AT(b,m+j-1,0))) {
-        case 0: { t = GLP_FR; break; } 
-        case 1: { t = GLP_LO; break; }
-        case 2: { t = GLP_UP; break; }
-        case 3: { t = GLP_DB; break; }
-       default: { t = GLP_FX; break; }  
-    }
-    glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
-    }
-    glp_term_out(0);
-    glp_simplex(lp, NULL);
-    sp[0] = glp_get_status(lp);
-    sp[1] = glp_get_obj_val(lp);
-    for (k=1; k<=n; k++) {
-        sp[k+1] = glp_get_col_prim(lp, k);
-    }
-    
-    glp_delete_prob(lp);
-    free(ia);
-    free(ja);
-    free(ar);
-
-    return 0;
-}
-
-/* ---------------------------------------------------- */
+/*-----------------------------------------------------*/
 
 int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
     glp_prob *lp;
@@ -108,25 +38,25 @@ int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
         //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
     }
     glp_load_matrix(lp, tot, ia, ja, ar);
-    
+
     int t;
     for (i=1;i<=m;i++) {
     switch((int)rint(AT(b,i-1,0))) {
-        case 0: { t = GLP_FR; break; } 
+        case 0: { t = GLP_FR; break; }
         case 1: { t = GLP_LO; break; }
         case 2: { t = GLP_UP; break; }
         case 3: { t = GLP_DB; break; }
-       default: { t = GLP_FX; break; }  
+       default: { t = GLP_FX; break; }
     }
     glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
     }
     for (j=1;j<=n;j++) {
     switch((int)rint(AT(b,m+j-1,0))) {
-        case 0: { t = GLP_FR; break; } 
+        case 0: { t = GLP_FR; break; }
         case 1: { t = GLP_LO; break; }
         case 2: { t = GLP_UP; break; }
         case 3: { t = GLP_DB; break; }
-       default: { t = GLP_FX; break; }  
+       default: { t = GLP_FX; break; }
     }
     glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
     }
@@ -144,4 +74,3 @@ int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
 
     return 0;
 }
-