diff --git a/cylp/cy/CyClpSimplex.pyx b/cylp/cy/CyClpSimplex.pyx
index c7d5cf2a..d2d6c0c9 100644
--- a/cylp/cy/CyClpSimplex.pyx
+++ b/cylp/cy/CyClpSimplex.pyx
@@ -150,7 +150,7 @@ cdef class CyClpSimplex:
                 # This does work in some versions of SciPy
                 # It would probably be OK if csr_matrixPlus didn't override
                 # __get_item__ to always cast the result back to csr_matrixPlus
-                # I'm not actually sure why the objective is stored as 
+                # I'm not actually sure why the objective is stored as
                 # csr_matrixPlus anyway... seems to not always be true.
                 #
                 #if isinstance(o, (sparse.coo_matrix,
@@ -292,6 +292,15 @@ cdef class CyClpSimplex:
         def __get__(self):
             return self.CppSelf.objectiveValue()
 
+    property primalColumnSolution:
+        '''
+        Raw Solution to the primal variables.
+
+        :rtype: Numpy array
+        '''
+        def __get__(self):
+            return <object>self.CppSelf.getPrimalColumnSolution()
+
     property primalVariableSolution:
         '''
         Solution to the primal variables.
@@ -388,6 +397,17 @@ cdef class CyClpSimplex:
                     ret = d
             return ret
 
+
+    property dualColumnSolution:
+        '''
+        reduced costs
+
+        :rtype: Numpy array
+        '''
+        def __get__(self):
+            return <object>self.CppSelf.getDualColumnSolution()
+
+
     property primalConstraintSolution:
         '''
         Slack variables' solution
diff --git a/cylp/tests/test_CyClpSimplex_CyLPModel.py b/cylp/tests/test_CyClpSimplex_CyLPModel.py
index 784229b9..7b47c149 100644
--- a/cylp/tests/test_CyClpSimplex_CyLPModel.py
+++ b/cylp/tests/test_CyClpSimplex_CyLPModel.py
@@ -262,7 +262,51 @@ def test_ArrayIndexing(self):
         self.assertTrue(abs(sol[1, 2, 3] - 1) <= 10**-6)
         self.assertTrue(abs(sol[1, 2, 5] - 1) <= 10**-6)
 
+    def test_primalAndDualColumnSolutions(self):
+        self.model = CyLPModel()
+        model = self.model
+
+        x = model.addVariable('x', 3)
+        y = model.addVariable('y', 2)
+
+        A = np.matrix([[1., 2., 0],[1., 0, 1.]])
+        B = np.matrix([[1., 0, 0], [0, 0, 1.]])
+        D = np.matrix([[1., 2.],[0, 1]])
+        a = CyLPArray([5, 2.5])
+        b = CyLPArray([4.2, 3])
+        x_u= CyLPArray([2., 3.5])
+
+        model.addConstraint(A*x <= a)
+        model.addConstraint(2 <= B * x + D * y <= b)
+        model.addConstraint(y >= 0)
+        model.addConstraint(1.1 <= x[1:3] <= x_u)
+
+        c = CyLPArray([1., -2., 3.])
+        model.objective = c * x + 2 * y.sum()
+
+        self.s = CyClpSimplex(model)
+        s = self.s
+
+        s.initialDualSolve()
+        self.assertEqual(s.getStatusCode(), 0)
+        self.assertEqual(round(s.objectiveValue, 4), 1.3000)
+
+        primal_sol = s.primalVariableSolution
+        dual_sol = s.dualVariableSolution
+
+        # assert that primalColumnSolution and primalVariableSolution are equivalent
+        # i.e. primalColumnSolution is single array containing all the primal variable solutions
+        self.assertTrue((s.primalColumnSolution == np.concatenate([primal_sol['x'], primal_sol['y']])).all())
+
+        # assert that dualColumnSolution and dualVariableSolution are equivalent
+        # i.e. dualColumnSolution is single array containing all the dual variable solutions
+        self.assertTrue((s.dualColumnSolution == np.concatenate([dual_sol['x'], dual_sol['y']])).all())
+
+        self.assertTrue((abs(primal_sol['x'] - np.array([0.2, 2.0 , 1.1]) ) <= 10**-6).all())
+        self.assertTrue((abs(primal_sol['y'] - np.array([0.0, 0.9]) ) <= 10**-6).all())
+
+        self.assertTrue((abs(dual_sol['x'] - np.array([0.0, -2.0,  3.0]) ) <= 10**-6).all())
+        self.assertTrue((abs(dual_sol['y'] - np.array([1.0, 0.0]) ) <= 10**-6).all())
 
 if __name__ == '__main__':
     unittest.main()
-