summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/glpk/lib/Numeric/LinearProgramming.hs55
-rw-r--r--packages/glpk/lib/Numeric/LinearProgramming/glpk.c83
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:
28 28
29prob = Maximize [4, -3, 2] 29prob = Maximize [4, -3, 2]
30 30
31constr1 = Sparse [ [2#1, 1#2] :<: 10 31constr1 = Sparse [ [2\#1, 1\#2] :<: 10
32 , [1#2, 5#3] :<: 20 32 , [1\#2, 5\#3] :<: 20
33 ] 33 ]
34 34
35\> simplex prob constr1 [] 35\> simplex prob constr1 []
@@ -52,6 +52,8 @@ Unbounded@
52 52
53The given bound for a variable completely replaces the default, 53The given bound for a variable completely replaces the default,
54so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@. 54so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@.
55Multiple bounds for a variable are not allowed, instead of
56@[i :>: a, i:<: b]@ use @i :&: (a,b)@.
55 57
56-} 58-}
57 59
@@ -65,11 +67,11 @@ module Numeric.LinearProgramming(
65 Solution(..) 67 Solution(..)
66) where 68) where
67 69
68import Numeric.LinearAlgebra 70import Numeric.LinearAlgebra hiding (i)
69import Data.Packed.Development 71import Data.Packed.Development
70import Foreign(Ptr,unsafePerformIO) 72import Foreign(Ptr,unsafePerformIO)
71import Foreign.C.Types(CInt) 73import Foreign.C.Types(CInt)
72import Data.List((\\),sortBy) 74import Data.List((\\),sortBy,nub)
73import Data.Function(on) 75import Data.Function(on)
74 76
75--import Debug.Trace 77--import Debug.Trace
@@ -82,7 +84,7 @@ import Data.Function(on)
82infixl 5 # 84infixl 5 #
83(#) = (,) 85(#) = (,)
84 86
85data Bound x = x :<: Double 87data Bound x = x :<: Double
86 | x :>: Double 88 | x :>: Double
87 | x :&: (Double,Double) 89 | x :&: (Double,Double)
88 | x :==: Double 90 | x :==: Double
@@ -107,11 +109,13 @@ type Bounds = [Bound Int]
107 109
108simplex :: Optimization -> Constraints -> Bounds -> Solution 110simplex :: Optimization -> Constraints -> Bounds -> Solution
109 111
110simplex opt (Dense []) bnds = simplex opt (Sparse [Free [0#1]]) bnds 112simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds
111simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds 113simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
112 114
113simplex opt (Dense constr) bnds = extract sg sol where 115simplex opt (Dense constr) bnds = extract sg sol where
114 sol = simplexDense (mkConstrD sz objfun constr) (mkBounds sz constr bnds) 116 sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
117 n = length objfun
118 m = length constr
115 (sz, sg, objfun) = adapt opt 119 (sz, sg, objfun) = adapt opt
116 120
117simplex opt (Sparse constr) bnds = extract sg sol where 121simplex opt (Sparse constr) bnds = extract sg sol where
@@ -178,46 +182,37 @@ mkBound2 b = (obj b, mkBound1 b)
178 182
179mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double 183mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double
180mkBounds n b1 b2 = fromLists (cb++vb) where 184mkBounds n b1 b2 = fromLists (cb++vb) where
181 gv = map obj b2 185 gv' = map obj b2
186 gv | nub gv' == gv' = gv'
187 | otherwise = error $ "simplex: duplicate bounds for vars " ++ show (gv'\\nub gv')
182 rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv 188 rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv
183 | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n 189 | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n
184 vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2 190 vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2
185 cb = map mkBound1 b1 191 cb = map mkBound1 b1
186 192
187mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double 193mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double
188mkConstrD n f b1 | ok = fromLists (f : cs) 194mkConstrD n f b1 | ok = fromLists (ob ++ co)
189 | otherwise = error $ "simplex: dense constraints require "++show n 195 | otherwise = error $ "simplex: dense constraints require "++show n
190 ++" variables, given " ++ show ls 196 ++" variables, given " ++ show ls
191 where 197 where
192 cs | null b1 = error $ "simplex: dense: empty constraints" 198 cs = map obj b1
193 | otherwise = map obj b1
194 ls = map length cs 199 ls = map length cs
195 ok = all (==n) ls 200 ok = all (==n) ls
201 den = fromLists cs
202 ob = map (([0,0]++).return) f
203 co = [[fromIntegral i, fromIntegral j,den@@>(i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]
196 204
197mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double 205mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double
198mkConstrS n objfun b1 = fromLists (ob ++ co) where 206mkConstrS n objfun b1 = fromLists (ob ++ co) where
199 ob = map (([0,0]++).return) objfun 207 ob = map (([0,0]++).return) objfun
200 co = concat $ zipWith f [1::Int ..] cs 208 co = concat $ zipWith f [1::Int ..] cs
201 cs | null b1 = error $ "simplex: sparse: empty constraints" 209 cs = map obj b1
202 | otherwise = map obj b1
203 f k = map (g k) 210 f k = map (g k)
204 g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c] 211 g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c]
205 | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n 212 | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n
206 213
207----------------------------------------------------- 214-----------------------------------------------------
208 215
209foreign import ccall "c_simplex_dense" c_simplex_dense
210 :: CInt -> CInt -> Ptr Double -- coeffs
211 -> CInt -> CInt -> Ptr Double -- bounds
212 -> CInt -> Ptr Double -- result
213 -> IO CInt -- exit code
214
215simplexDense :: Matrix Double -> Matrix Double -> Vector Double
216simplexDense c b = unsafePerformIO $ do
217 s <- createVector (2+cols c)
218 app3 c_simplex_dense mat (cmat c) mat (cmat b) vec s "c_simplex_dense"
219 return s
220
221foreign import ccall "c_simplex_sparse" c_simplex_sparse 216foreign import ccall "c_simplex_sparse" c_simplex_sparse
222 :: CInt -> CInt -- rows and cols 217 :: CInt -> CInt -- rows and cols
223 -> CInt -> CInt -> Ptr Double -- coeffs 218 -> CInt -> CInt -> Ptr Double -- coeffs
@@ -241,14 +236,6 @@ glpFX = 4
241 236
242{- Raw format of coeffs 237{- Raw format of coeffs
243 238
244simplexDense
245
246((3+1) >< 3)
247 [ 10, 6, 4
248 , 1, 1, 1
249 , 10, 4, 5
250 , 2, 2, 6 :: Double]
251
252simplexSparse 239simplexSparse
253 240
254(12><3) 241(12><3)
@@ -264,7 +251,7 @@ simplexSparse
264 , 3.0, 1.0, 2.0 251 , 3.0, 1.0, 2.0
265 , 3.0, 2.0, 2.0 252 , 3.0, 2.0, 2.0
266 , 3.0, 3.0, 6.0 ] 253 , 3.0, 3.0, 6.0 ]
267 254
268bounds = (6><3) 255bounds = (6><3)
269 [ glpUP,0,100 256 [ glpUP,0,100
270 , glpUP,0,600 257 , 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 @@
1typedef struct { double r, i; } doublecomplex;
2
3#define DVEC(A) int A##n, double*A##p 1#define DVEC(A) int A##n, double*A##p
4#define CVEC(A) int A##n, doublecomplex*A##p
5#define DMAT(A) int A##r, int A##c, double*A##p 2#define DMAT(A) int A##r, int A##c, double*A##p
6#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
7 3
8#define AT(M,r,co) (M##p[(r)*M##c+(co)]) 4#define AT(M,r,co) (M##p[(r)*M##c+(co)])
9 5
10/*-----------------------------------------------------*/
11
12#include <stdlib.h> 6#include <stdlib.h>
13#include <stdio.h> 7#include <stdio.h>
14#include <glpk.h> 8#include <glpk.h>
15#include <math.h> 9#include <math.h>
16 10
17int c_simplex_dense(DMAT(c), DMAT(b), DVEC(s)) { 11/*-----------------------------------------------------*/
18 glp_prob *lp;
19 lp = glp_create_prob();
20 glp_set_obj_dir(lp, GLP_MAX);
21
22 int m = cr-1;
23 int n = cc;
24 glp_add_rows(lp, m);
25 glp_add_cols(lp, n);
26 int * ia = malloc((1+m*n)*sizeof(int));
27 int * ja = malloc((1+m*n)*sizeof(int));
28 double * ar = malloc((1+m*n)*sizeof(double));
29 int i,j,k;
30 k=0;
31 for (i=1;i<=m;i++) {
32 for(j=1;j<=n;j++) {
33 k++;
34 ia[k] = i;
35 ja[k] = j;
36 ar[k] = AT(c,i,j-1);
37 //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
38 }
39 }
40 glp_load_matrix(lp, k, ia, ja, ar);
41 for (j=1;j<=n;j++) {
42 glp_set_obj_coef(lp, j, AT(c,0,j-1));
43 }
44 int t;
45 for (i=1;i<=m;i++) {
46 switch((int)rint(AT(b,i-1,0))) {
47 case 0: { t = GLP_FR; break; }
48 case 1: { t = GLP_LO; break; }
49 case 2: { t = GLP_UP; break; }
50 case 3: { t = GLP_DB; break; }
51 default: { t = GLP_FX; break; }
52 }
53 glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
54 }
55 for (j=1;j<=n;j++) {
56 switch((int)rint(AT(b,m+j-1,0))) {
57 case 0: { t = GLP_FR; break; }
58 case 1: { t = GLP_LO; break; }
59 case 2: { t = GLP_UP; break; }
60 case 3: { t = GLP_DB; break; }
61 default: { t = GLP_FX; break; }
62 }
63 glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
64 }
65 glp_term_out(0);
66 glp_simplex(lp, NULL);
67 sp[0] = glp_get_status(lp);
68 sp[1] = glp_get_obj_val(lp);
69 for (k=1; k<=n; k++) {
70 sp[k+1] = glp_get_col_prim(lp, k);
71 }
72
73 glp_delete_prob(lp);
74 free(ia);
75 free(ja);
76 free(ar);
77
78 return 0;
79}
80
81/* ---------------------------------------------------- */
82 12
83int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) { 13int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
84 glp_prob *lp; 14 glp_prob *lp;
@@ -108,25 +38,25 @@ int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
108 //printf("%d %d %f\n",ia[k],ja[k],ar[k]); 38 //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
109 } 39 }
110 glp_load_matrix(lp, tot, ia, ja, ar); 40 glp_load_matrix(lp, tot, ia, ja, ar);
111 41
112 int t; 42 int t;
113 for (i=1;i<=m;i++) { 43 for (i=1;i<=m;i++) {
114 switch((int)rint(AT(b,i-1,0))) { 44 switch((int)rint(AT(b,i-1,0))) {
115 case 0: { t = GLP_FR; break; } 45 case 0: { t = GLP_FR; break; }
116 case 1: { t = GLP_LO; break; } 46 case 1: { t = GLP_LO; break; }
117 case 2: { t = GLP_UP; break; } 47 case 2: { t = GLP_UP; break; }
118 case 3: { t = GLP_DB; break; } 48 case 3: { t = GLP_DB; break; }
119 default: { t = GLP_FX; break; } 49 default: { t = GLP_FX; break; }
120 } 50 }
121 glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2)); 51 glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
122 } 52 }
123 for (j=1;j<=n;j++) { 53 for (j=1;j<=n;j++) {
124 switch((int)rint(AT(b,m+j-1,0))) { 54 switch((int)rint(AT(b,m+j-1,0))) {
125 case 0: { t = GLP_FR; break; } 55 case 0: { t = GLP_FR; break; }
126 case 1: { t = GLP_LO; break; } 56 case 1: { t = GLP_LO; break; }
127 case 2: { t = GLP_UP; break; } 57 case 2: { t = GLP_UP; break; }
128 case 3: { t = GLP_DB; break; } 58 case 3: { t = GLP_DB; break; }
129 default: { t = GLP_FX; break; } 59 default: { t = GLP_FX; break; }
130 } 60 }
131 glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2)); 61 glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
132 } 62 }
@@ -144,4 +74,3 @@ int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
144 74
145 return 0; 75 return 0;
146} 76}
147