diff options
author | Alberto Ruiz <aruiz@um.es> | 2010-02-24 09:03:40 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2010-02-24 09:03:40 +0000 |
commit | 4ebde988645cb1d2d60f8df8605ec8254713a9d9 (patch) | |
tree | 3fd001a92921dcc553e01c8c32849f5e8135247a | |
parent | e0605467b60f65478d0f5cc8f82ba14a99168f7b (diff) |
removed simplexDense
-rw-r--r-- | packages/glpk/lib/Numeric/LinearProgramming.hs | 55 | ||||
-rw-r--r-- | packages/glpk/lib/Numeric/LinearProgramming/glpk.c | 83 |
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 | ||
29 | prob = Maximize [4, -3, 2] | 29 | prob = Maximize [4, -3, 2] |
30 | 30 | ||
31 | constr1 = Sparse [ [2#1, 1#2] :<: 10 | 31 | constr1 = 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 | ||
53 | The given bound for a variable completely replaces the default, | 53 | The given bound for a variable completely replaces the default, |
54 | so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@. | 54 | so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@. |
55 | Multiple 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 | ||
68 | import Numeric.LinearAlgebra | 70 | import Numeric.LinearAlgebra hiding (i) |
69 | import Data.Packed.Development | 71 | import Data.Packed.Development |
70 | import Foreign(Ptr,unsafePerformIO) | 72 | import Foreign(Ptr,unsafePerformIO) |
71 | import Foreign.C.Types(CInt) | 73 | import Foreign.C.Types(CInt) |
72 | import Data.List((\\),sortBy) | 74 | import Data.List((\\),sortBy,nub) |
73 | import Data.Function(on) | 75 | import Data.Function(on) |
74 | 76 | ||
75 | --import Debug.Trace | 77 | --import Debug.Trace |
@@ -82,7 +84,7 @@ import Data.Function(on) | |||
82 | infixl 5 # | 84 | infixl 5 # |
83 | (#) = (,) | 85 | (#) = (,) |
84 | 86 | ||
85 | data Bound x = x :<: Double | 87 | data 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 | ||
108 | simplex :: Optimization -> Constraints -> Bounds -> Solution | 110 | simplex :: Optimization -> Constraints -> Bounds -> Solution |
109 | 111 | ||
110 | simplex opt (Dense []) bnds = simplex opt (Sparse [Free [0#1]]) bnds | 112 | simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds |
111 | simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds | 113 | simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds |
112 | 114 | ||
113 | simplex opt (Dense constr) bnds = extract sg sol where | 115 | simplex 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 | ||
117 | simplex opt (Sparse constr) bnds = extract sg sol where | 121 | simplex opt (Sparse constr) bnds = extract sg sol where |
@@ -178,46 +182,37 @@ mkBound2 b = (obj b, mkBound1 b) | |||
178 | 182 | ||
179 | mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double | 183 | mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double |
180 | mkBounds n b1 b2 = fromLists (cb++vb) where | 184 | mkBounds 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 | ||
187 | mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double | 193 | mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double |
188 | mkConstrD n f b1 | ok = fromLists (f : cs) | 194 | mkConstrD 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 | ||
197 | mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double | 205 | mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double |
198 | mkConstrS n objfun b1 = fromLists (ob ++ co) where | 206 | mkConstrS 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 | ||
209 | foreign 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 | |||
215 | simplexDense :: Matrix Double -> Matrix Double -> Vector Double | ||
216 | simplexDense 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 | |||
221 | foreign import ccall "c_simplex_sparse" c_simplex_sparse | 216 | foreign 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 | ||
244 | simplexDense | ||
245 | |||
246 | ((3+1) >< 3) | ||
247 | [ 10, 6, 4 | ||
248 | , 1, 1, 1 | ||
249 | , 10, 4, 5 | ||
250 | , 2, 2, 6 :: Double] | ||
251 | |||
252 | simplexSparse | 239 | simplexSparse |
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 | ||
268 | bounds = (6><3) | 255 | bounds = (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 @@ | |||
1 | typedef 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 | ||
17 | int 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 | ||
83 | int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) { | 13 | int 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 | |||