From 4ebde988645cb1d2d60f8df8605ec8254713a9d9 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 24 Feb 2010 09:03:40 +0000 Subject: removed simplexDense --- packages/glpk/lib/Numeric/LinearProgramming.hs | 55 ++++++-------- packages/glpk/lib/Numeric/LinearProgramming/glpk.c | 83 ++-------------------- 2 files changed, 27 insertions(+), 111 deletions(-) (limited to 'packages') 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 #include #include #include -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; } - -- cgit v1.2.3