diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -97,3 +97,7 @@
 
 0.4.1.1:
 		remove Config.hs and configure
+
+0.4.1.2:
+		update for hmatrix 0.17
+
diff --git a/hmatrix-gsl-stats.cabal b/hmatrix-gsl-stats.cabal
--- a/hmatrix-gsl-stats.cabal
+++ b/hmatrix-gsl-stats.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-gsl-stats
-Version:            0.4.1.1
+Version:            0.4.1.2
 License:            BSD3
 License-file:       LICENSE
 Copyright:          (c) A.V.H. McPhail 2010, 2011, 2013, 2015
@@ -17,7 +17,7 @@
      .
      The vector type is Data.Vector.Storable from the 'vector' package.
 Category:           Math, Statistics, FFI
-tested-with:        GHC ==7.10.1
+tested-with:        GHC ==7.10.2
 
 cabal-version:      >=1.10.1.0
 
@@ -35,8 +35,9 @@
     Default-language:   Haskell2010
 
     Build-Depends:      base >= 4 && < 5, binary,
+                        vector >= 0.11,
                         storable-complex,
-                        hmatrix >= 0.16
+                        hmatrix >= 0.17
 			
     Default-Extensions: ForeignFunctionInterface
 
diff --git a/lib/Numeric/GSL/Distribution/Continuous.hs b/lib/Numeric/GSL/Distribution/Continuous.hs
--- a/lib/Numeric/GSL/Distribution/Continuous.hs
+++ b/lib/Numeric/GSL/Distribution/Continuous.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE TypeFamilies #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Distribution.Continuous
@@ -32,12 +33,9 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
---import Data.Packed.Matrix hiding(toLists)
-import Data.Packed.Development
---import Data.Packed.Development(createVector,vec,app1,app2,app3,app4)
 
---import Numeric.LinearAlgebra.Linear
+import Numeric.LinearAlgebra.Data
+import Numeric.LinearAlgebra.Devel
 
 --import Control.Monad(when)
 
@@ -62,6 +60,12 @@
 
 -----------------------------------------------------------------------------
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 data ZeroParamDist = Landau
                    deriving Enum
 
@@ -124,7 +128,7 @@
             -> Vector Double   -- ^ result
 random_0p_v  d s l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_random_zero_param_v (fromIntegral s) (fromei d)) vec r "random_0p_v"
+   (distribution_random_zero_param_v (fromIntegral s) (fromei d)) # r #| "random_0p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -175,7 +179,7 @@
             -> Vector Double   -- ^ result
 random_1p_v d s p l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_random_one_param_v (fromIntegral s) (fromei d) p) vec r "random_1p_v"
+   (distribution_random_one_param_v (fromIntegral s) (fromei d) p) # r #| "random_1p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -225,7 +229,7 @@
             -> Vector Double   -- ^ result
 random_2p_v d s p1 p2 l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_random_two_param_v (fromIntegral s) (fromei d) p1 p2) vec r "random_2p_v"
+   (distribution_random_two_param_v (fromIntegral s) (fromei d) p1 p2) # r #| "random_2p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -291,7 +295,7 @@
             -> Vector Double   -- ^ result
 random_3p_v d s p1 p2 p3 l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_random_three_param_v (fromIntegral s) (fromei d) p1 p2 p3) vec r "random_3p_v"
+   (distribution_random_three_param_v (fromIntegral s) (fromei d) p1 p2 p3) # r #| "random_3p_v"
    return r
 
 -- | probability of a variate take a value outside the argument
@@ -320,8 +324,8 @@
           -> Vector Double   -- ^ parameters
           -> Vector Double   -- ^ result
 random_mp d s p = unsafePerformIO $ do
-                  r <- createVector $ dim p
-                  app2 (distribution_random_multi_param (fromIntegral s) (fromei d)) vec p vec r "random_mp"
+                  r <- createVector $ size p
+                  (distribution_random_multi_param (fromIntegral s) (fromei d)) # p # r #| "random_mp"
                   return r
 
 -- | draw a sample from a multi parameter distribution
@@ -330,9 +334,9 @@
             -> Vector Double      -- ^ parameters
             -> IO (Vector Double) -- ^ result
 random_mp_s (RNG rng) d p = do
-  let l = dim p
+  let l = size p
   r <- createVector l
-  withForeignPtr rng $ \rg -> app2 (distribution_random_multi_param_s rg (fromei d)) vec p vec r "random_mp_s"
+  withForeignPtr rng $ \rg -> (distribution_random_multi_param_s rg (fromei d)) # p # r #| "random_mp_s"
   return r
   
 -- | probability of a variate take a value outside the argument
@@ -347,7 +351,7 @@
     where density_only f' d' p' q' = if f' /= Density
                                               then error "distribution has no CDF"
                                               else alloca $ \r -> do
-                                                                  app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p' vec q' "density_mp"
+                                                                  (distribution_dist_multi_param (fromei f') (fromei d') r) # p' # q' #| "density_mp"
                                                                   r' <- peek r
                                                                   return r'
 
@@ -398,7 +402,7 @@
 random_biv_v d s p1 p2 p3 l = unsafePerformIO $ do
    r1 <- createVector l
    r2 <- createVector l
-   app2 (distribution_random_bivariate_v (fromIntegral s) (fromei d) p1 p2 p3) vec r1 vec r2 "random_biv_v"
+   (distribution_random_bivariate_v (fromIntegral s) (fromei d) p1 p2 p3) # r1 # r2 #| "random_biv_v"
    return (r1,r2)
 
 -- | probability of a variate take a value outside the argument
@@ -429,7 +433,7 @@
                  -> Vector Double -- result
 spherical_vector s vs = unsafePerformIO $ do
                         r <- createVector vs
-                        app1 (distribution_spherical_vector (fromIntegral s)) vec r "spherical_vector"
+                        (distribution_spherical_vector (fromIntegral s)) # r #| "spherical_vector"
                         return r
 
 foreign import ccall "distribution-aux.h spherical_vector" distribution_spherical_vector :: CInt -> CInt -> Ptr Double -> IO CInt
diff --git a/lib/Numeric/GSL/Distribution/Discrete.hs b/lib/Numeric/GSL/Distribution/Discrete.hs
--- a/lib/Numeric/GSL/Distribution/Discrete.hs
+++ b/lib/Numeric/GSL/Distribution/Discrete.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE TypeFamilies #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Distribution.Discrete
@@ -27,12 +28,10 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
---import Data.Packed.Matrix hiding(toLists)
-import Data.Packed.Development
-import Data.Word
+import Numeric.LinearAlgebra.Data
+import Numeric.LinearAlgebra.Devel
 
---import Numeric.LinearAlgebra.Linear
+import Data.Vector.Storable
 
 --import Control.Monad(when)
 
@@ -48,7 +47,7 @@
 --import GHC.Base
 --import GHC.IOBase
 
---import Prelude hiding(reverse)
+import Prelude hiding(map)
 
 import Numeric.GSL.Distribution.Common
 import Numeric.GSL.Distribution.Internal
@@ -57,6 +56,12 @@
 
 -----------------------------------------------------------------------------
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 data OneParamDist = Poisson     -- ^ mean
                   | Bernoulli   -- ^ probability
                   | Geometric   -- ^ probability
@@ -111,8 +116,8 @@
             -> Vector Word32   -- ^ result
 random_1p_v d s p l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_discrete_one_param_v (fromIntegral s) (fromei d) p) vec r "random_1p_v"
-   return (mapVector fromIntegral r)
+   (distribution_discrete_one_param_v (fromIntegral s) (fromei d) p) # r #| "random_1p_v"
+   return (map fromIntegral r)
 
 -- | probability of a variate take a value outside the argument
 density_1p :: OneParamDist   -- ^ density type
@@ -174,8 +179,8 @@
             -> Vector Word32   -- ^ result
 random_2p_v d s p1 p2 l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_discrete_two_param_v (fromIntegral s) (fromei d) p1 (fromIntegral p2)) vec r "random_2p_v"
-   return (mapVector fromIntegral r)
+   (distribution_discrete_two_param_v (fromIntegral s) (fromei d) p1 (fromIntegral p2)) # r #| "random_2p_v"
+   return (map fromIntegral r)
 
 -- | probability of a variate take a value outside the argument
 density_2p :: TwoParamDist   -- ^ density type
@@ -234,8 +239,8 @@
             -> Vector Word32   -- ^ result
 random_3p_v d s p1 p2 p3 l = unsafePerformIO $ do
    r <- createVector l
-   app1 (distribution_discrete_three_param_v (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3)) vec r "random_3p_v"
-   return (mapVector fromIntegral r)
+   (distribution_discrete_three_param_v (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3)) # r #| "random_3p_v"
+   return (map fromIntegral r)
 
 -- | probability of a variate take a value outside the argument
 density_3p :: ThreeParamDist   -- ^ density type
@@ -266,9 +271,9 @@
           -> Vector Double   -- ^ parameters
           -> Vector Word32   -- ^ result
 random_mp d s t p = unsafePerformIO $ do
-                    r <- createVector $ dim p
-                    app2 (distribution_discrete_multi_param (fromIntegral s) (fromei d) (fromIntegral t)) vec p vec r "random_mp"
-                    return $ mapVector (\x -> (fromIntegral x) :: Word32) r
+                    r <- createVector $ size p
+                    (distribution_discrete_multi_param (fromIntegral s) (fromei d) (fromIntegral t)) # p # r #| "random_mp"
+                    return $ map (\x -> (fromIntegral x) :: Word32) r
 
 -- | draw a sample from a multi-parameter distribution
 random_mp_s :: RNG           -- ^ RNG
@@ -277,9 +282,9 @@
             -> Vector Double   -- ^ parameters
             -> IO (Vector Word32) -- ^ result
 random_mp_s (RNG rng) d t p = withForeignPtr rng $ \rg -> do
-  r <- createVector $ dim p
-  app2 (distribution_discrete_multi_param_s rg (fromei d) (fromIntegral t)) vec p vec r "random_mp_s"
-  return $ mapVector (\x -> (fromIntegral x) :: Word32) r
+  r <- createVector $ size p
+  (distribution_discrete_multi_param_s rg (fromei d) (fromIntegral t)) # p # r #| "random_mp_s"
+  return $ map (\x -> (fromIntegral x) :: Word32) r
 
 -- | probability of a variate take a value outside the argument
 density_mp :: MultiParamDist   -- ^ density type
@@ -293,7 +298,7 @@
     where density_only f' d' p' q' = if f' /= Density
                                               then error "distribution has no CDF"
                                               else alloca $ \r -> do
-                                                                  app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p' vec (mapVector (\x -> (fromIntegral x) :: CUInt) q') "density_mp"
+                                                                  (distribution_dist_multi_param (fromei f') (fromei d') r) # p' # (map (\x -> (fromIntegral x) :: CUInt) q') #| "density_mp"
                                                                   r' <- peek r
                                                                   return r'
 
diff --git a/lib/Numeric/GSL/Fitting/Linear.hs b/lib/Numeric/GSL/Fitting/Linear.hs
--- a/lib/Numeric/GSL/Fitting/Linear.hs
+++ b/lib/Numeric/GSL/Fitting/Linear.hs
@@ -21,9 +21,8 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
-import Data.Packed.Matrix
-import Data.Packed.Development
+import Numeric.LinearAlgebra.Data
+import Numeric.LinearAlgebra.Devel
 
 --import Numeric.LinearAlgebra.Linear
 
@@ -46,6 +45,12 @@
 
 -----------------------------------------------------------------------------
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 -- | fits the model Y = C X
 linear :: Vector Double    -- ^ x data
        -> Vector Double    -- ^ y data
@@ -57,7 +62,7 @@
                        alloca $ \cov00 -> 
                            alloca $ \cov01 -> 
                                alloca $ \cov11 -> do
-                                                  app2 (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) vec x vec y "linear"
+                                                  (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) # x # y #| "linear"
                                                   c0' <- peek c0
                                                   c1' <- peek c1
                                                   cov00' <- peek cov00
@@ -85,7 +90,7 @@
                        alloca $ \cov00 -> 
                            alloca $ \cov01 -> 
                                alloca $ \cov11 -> do
-                                                  app3 (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) vec x vec w vec y "linear_w"
+                                                  (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) # x # w # y #| "linear_w"
                                                   c0' <- peek c0
                                                   c1' <- peek c1
                                                   cov00' <- peek cov00
@@ -134,7 +139,7 @@
                cov <- createMatrix RowMajor p p
                c <- createVector p
                alloca$ \chi_sq -> do
-                   app4 (fitting_multifit chi_sq) mat (cmat x) vec y vec c mat cov "multifit"
+                   (fitting_multifit chi_sq) # (cmat x) # y # c # cov #| "multifit"
                    chi_sq' <- peek chi_sq
                    return (c,cov,chi_sq')
 
@@ -154,7 +159,7 @@
                    cov <- createMatrix RowMajor p p
                    c <- createVector p
                    alloca$ \chi_sq -> do
-                       app5 (fitting_multifit_w chi_sq) mat (cmat x) vec w vec y vec c mat cov "multifit"
+                       (fitting_multifit_w chi_sq) # (cmat x) # w # y # c # cov #| "multifit"
                        chi_sq' <- peek chi_sq
                        return (c,cov,chi_sq')
 
@@ -172,7 +177,7 @@
 multifit_est x c cov = unsafePerformIO $ do
                        alloca $ \y ->
                            alloca $ \e -> do
-                               app3 (fitting_multifit_est y e) vec x vec c mat cov "multifit_estimate"
+                               (fitting_multifit_est y e) # x # c # cov #| "multifit_estimate"
                                y' <- peek y
                                e' <- peek e
                                return (y',e')
diff --git a/lib/Numeric/GSL/Histogram.hs b/lib/Numeric/GSL/Histogram.hs
--- a/lib/Numeric/GSL/Histogram.hs
+++ b/lib/Numeric/GSL/Histogram.hs
@@ -47,12 +47,11 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
---import Data.Packed.Matrix hiding(toLists)
-import Data.Packed.Development
+import Numeric.LinearAlgebra.Data hiding(find)
+import Numeric.LinearAlgebra.Devel
 
 --import Numeric.LinearAlgebra.Algorithms hiding (multiply)
-import Numeric.LinearAlgebra hiding (multiply,add,divide,scale,find)
+--import Numeric.LinearAlgebra hiding (multiply,add,divide,scale,find)
 --import Numeric.Container 
 
 --import Control.Monad
@@ -78,6 +77,12 @@
 -----------------------------------------------------------------------------
 
 data Hist
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 type HistHandle = Ptr Hist
 -- | A histogram structure
 data Histogram = H { hdim :: {-# UNPACK #-} !Int -- ^ number of bins
@@ -128,10 +133,10 @@
 -- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1)
 fromRangesIO :: Vector Double -> IO Histogram
 fromRangesIO v = do
-               let sz = fromIntegral $ dim v - 1
+               let sz = fromIntegral $ size v - 1
                h <- histogram_new sz
                h' <- newForeignPtr histogram_free h
-               app1 (\d p -> withForeignPtr h' $ \f -> histogram_set_ranges f (fromIntegral d) p) vec v "fromRanges"
+               (\d p -> withForeignPtr h' $ \f -> histogram_set_ranges f (fromIntegral d) p) # v #| "fromRanges"
                return $ H (fromIntegral sz) h'
 
 -- | create a histogram with n bins and lower and upper limits
@@ -193,7 +198,7 @@
 fromVectors :: Vector Double -> Vector Double -> Histogram
 fromVectors r b = unsafePerformIO $ do
                   h@(H _ h') <- fromRangesIO r
-                  app2 (\rs r' bs b' -> withForeignPtr h' $ \h'' -> histogram_from_vectors h'' rs r' bs b') vec r vec b "fromVectors"
+                  (\rs r' bs b' -> withForeignPtr h' $ \h'' -> histogram_from_vectors h'' rs r' bs b') # r # b #| "fromVectors"
                   return h
                   
 foreign import ccall "gsl-histogram.h from_vectors" histogram_from_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
@@ -203,7 +208,7 @@
 toVectors (H b h) = unsafePerformIO $ do
                     rs <- createVector (b+1)
                     bs <- createVector b
-                    app2 (\s1 p1 s2 p2 -> withForeignPtr h $ \f -> histogram_to_vectors f s1 p1 s2 p2) vec rs vec bs "toVectors"
+                    (\s1 p1 s2 p2 -> withForeignPtr h $ \f -> histogram_to_vectors f s1 p1 s2 p2) # rs # bs #| "toVectors"
                     return (rs,bs)
 
 foreign import ccall "gsl-histogram.h to_vectors" histogram_to_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
@@ -260,7 +265,7 @@
 -- | add 1.0 to the correct bin for each element of the vector, fails silently if the value is outside the range
 incrementVectorIO :: Histogram -> Vector Double -> IO ()
 incrementVectorIO (H _ h) v = do
-                             app1 (\s p -> withForeignPtr h (\f -> histogram_increment_vector f s p)) vec v "incrementVector"
+                             (\s p -> withForeignPtr h (\f -> histogram_increment_vector f s p)) # v #| "incrementVector"
                              return ()
 
 -- | adds the weight (second Double) to the bin appropriate for the value (first Double)
@@ -276,7 +281,7 @@
 -- | add the weight (second vector) to the correct bin for each element of the first vector, fails silently if the value is outside the range
 accumulateVectorIO :: Histogram -> Vector Double -> Vector Double -> IO ()
 accumulateVectorIO (H _ h) v w = do
-                                app2 (\s1 p1 s2 p2 -> withForeignPtr h (\f -> histogram_accumulate_vector f s1 p1 s2 p2)) vec v vec w "accumulateVector"
+                                (\s1 p1 s2 p2 -> withForeignPtr h (\f -> histogram_accumulate_vector f s1 p1 s2 p2)) # v # w #| "accumulateVector"
                                 return ()
 
 -- | returns the contents of the i-th bin
@@ -347,8 +352,8 @@
 
 -- | find the number of occurences for each element of the input vector
 count :: Histogram -> Vector Double -> Vector Double
-count h = mapVector (\x -> let Just x' = find h x
-                           in getBin h x')
+count h = cmap (\x -> let Just x' = find h x
+                      in getBin h x')
 
 -- | find the probability of occurring for each element of the input vector
 prob :: Histogram -> Vector Double -> Vector Double
diff --git a/lib/Numeric/GSL/Histogram2D.hs b/lib/Numeric/GSL/Histogram2D.hs
--- a/lib/Numeric/GSL/Histogram2D.hs
+++ b/lib/Numeric/GSL/Histogram2D.hs
@@ -48,14 +48,14 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
-import Data.Packed.Matrix
-import Data.Packed.Development
-import Data.Packed()
+import Numeric.LinearAlgebra.Data hiding(find)
+import Numeric.LinearAlgebra.Devel
 
+import Data.Vector.Storable hiding(Vector,sum,find,mapM_)
+
 --import Numeric.LinearAlgebra.Algorithms hiding (multiply)
 --import Numeric.LinearAlgebra.Linear hiding (multiply,add,divide,scale)
-import Numeric.LinearAlgebra hiding (multiply,add,divide,scale,find)
+--import Numeric.LinearAlgebra hiding (multiply,add,divide,scale,find)
 --import Numeric.Container 
 
 --import Control.Monad
@@ -77,12 +77,18 @@
 --import GHC.Base
 --import GHC.IOBase
 
-import Prelude hiding(subtract,sum)
+import Prelude hiding(subtract,sum,length)
 
 import System.IO.Unsafe(unsafePerformIO)
 
 -----------------------------------------------------------------------------
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 instance (Storable a) => Storable (a,a) where
     sizeOf z        = 2 * sizeOf (fst z)
     alignment z     = alignment (fst z)
@@ -150,11 +156,11 @@
 -- | create a histogram with n bins from ranges (x0->x1),(x1->x2),..,(xn->xn+1) and y0,..,yn+1
 fromRangesIO :: Vector Double -> Vector Double -> IO Histogram2D
 fromRangesIO v w = do
-                   let sx = fromIntegral $ dim v - 1
-                   let sy = fromIntegral $ dim w - 1
+                   let sx = fromIntegral $ size v - 1
+                   let sy = fromIntegral $ size w - 1
                    h <- histogram2d_new sx sy
                    h' <- newForeignPtr histogram2d_free h
-                   app2 (\xs xp ys yp -> withForeignPtr h' (\f -> histogram2d_set_ranges f xp xs yp ys)) vec v vec w "fromRanges"
+                   (\xs xp ys yp -> withForeignPtr h' (\f -> histogram2d_set_ranges f xp xs yp ys)) # v # w #| "fromRanges"
                    return $ H (fromIntegral sx) (fromIntegral sy) h'
 
 -- | create a histogram with n bins and lower and upper limits
@@ -215,7 +221,7 @@
                     rx <- createVector (bx+1)
                     ry <- createVector (by+1)
                     bs <- createMatrix RowMajor bx by
-                    app3 (\s1 p1 s2 p2 sx sy p -> withForeignPtr h $ \f -> histogram_to_matrix f s1 p1 s2 p2 sx sy p) vec rx vec ry mat bs "toMatrix"
+                    (\s1 p1 s2 p2 sx sy p -> withForeignPtr h $ \f -> histogram_to_matrix f s1 p1 s2 p2 sx sy p) # rx # ry # bs #| "toMatrix"
                     return (rx,ry,bs)
 
 foreign import ccall "histogram-aux.h to_matrix" histogram_to_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
@@ -237,7 +243,7 @@
            -> Histogram2D              -- ^result
 fromMatrix x y w = unsafePerformIO $ do
                    h@(H _ _ h') <- fromRangesIO x y
-                   app3 (\xs x' ys y' rs cs b -> withForeignPtr h' $ \h'' -> histogram_from_matrix h'' xs x' ys y' rs cs b) vec x vec y mat (cmat w) "fromMatrix"
+                   (\xs x' ys y' rs cs b -> withForeignPtr h' $ \h'' -> histogram_from_matrix h'' xs x' ys y' rs cs b) # x # y # (cmat w) #| "fromMatrix"
                    return h
 
 foreign import ccall "histogram-aux.h from_matrix" histogram_from_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
@@ -295,7 +301,7 @@
 -- | add 1.0 to the correct bin for each element of the vector pair, fails silently if the value is outside the range
 incrementVectorIO :: Histogram2D -> Vector Double -> Vector Double -> IO ()
 incrementVectorIO (H _ _ h) x y = do
-                                app2 (\xs xp ys yp -> withForeignPtr h (\f -> histogram2d_increment_matrix f xs xp ys yp)) vec x vec y "incrementVector"
+                                (\xs xp ys yp -> withForeignPtr h (\f -> histogram2d_increment_matrix f xs xp ys yp)) # x # y #| "incrementVector"
                                 return ()
 
 -- | add 1.0  to the correct bin for each element of the list, fails silently if the value is outside the range
@@ -311,7 +317,7 @@
 -- | add the weight (third) to the correct bin for each vector pair element, fails silently if the value is outside the range
 accumulateVectorIO :: Histogram2D -> Vector Double -> Vector Double -> Vector Double -> IO ()
 accumulateVectorIO (H _ _ h) x y w = do
-                                app3 (\xs xp ys yp ws wp -> withForeignPtr h (\f -> histogram2d_accumulate_matrix f xs xp ys yp ws wp)) vec x vec y vec w "accumulateVector"
+                                (\xs xp ys yp ws wp -> withForeignPtr h (\f -> histogram2d_accumulate_matrix f xs xp ys yp ws wp)) # x # y # w #| "accumulateVector"
                                 return ()
 
 -- | add the weight (snd) to the correct bin for each (fst) element of the list, fails silently if the value is outside the range
@@ -415,8 +421,8 @@
 -- | find the number of occurences for each element of the input vector
 count :: Histogram2D -> (Vector Double, Vector Double) -> Vector Double
 count (H _ _ h) (x,y) = unsafePerformIO $ do
-               r <- createVector $ dim x
-               app3 (\xs' x' ys' y' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count h' xs' x' ys' y' rs' r') vec x vec y vec r "histogram2d_count"
+               r <- createVector $ size x
+               (\xs' x' ys' y' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count h' xs' x' ys' y' rs' r') # x # y # r #| "histogram2d_count"
                return r
 
 foreign import ccall "histogram-aux.h hist2d_count" histogram2d_count :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
@@ -428,8 +434,8 @@
 -- | find the number of occurences for each element of the input vector
 countPaired :: Histogram2D -> Vector (Double,Double) -> Vector Double
 countPaired (H _ _ h) x = unsafePerformIO $ do
-               r <- createVector $ dim x
-               app2 (\xs' x' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count_pair h' xs' x' rs' r') vec x vec r "histogram2d_count_pair"
+               r <- createVector $ length x
+               (\xs' x' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count_pair h' xs' x' rs' r') # x # r #| "histogram2d_count_pair"
                return r
 
 foreign import ccall "histogram-aux.h hist2d_count_pair" histogram2d_count_pair :: Hist2DHandle -> CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> IO CInt
@@ -654,9 +660,9 @@
 {-
 unzipPair :: Vector (Double,Double) -> (Vector Double,Vector Double)
 unzipPair v = unsafePerformIO $ do
-              a <- createVector $ dim v
-              b <- createVector $ dim v
-              app3 histogram2d_unzip_double_pair vec v vec a vec b "unzipPair"
+              a <- createVector $ size v
+              b <- createVector $ size v
+              histogram2d_unzip_double_pair # v # a # b #| "unzipPair"
               return (a,b)
 
 foreign import ccall "histogram-aux.h unzip_double_pair" histogram2d_unzip_double_pair :: CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
diff --git a/lib/Numeric/GSL/Permutation.hs b/lib/Numeric/GSL/Permutation.hs
--- a/lib/Numeric/GSL/Permutation.hs
+++ b/lib/Numeric/GSL/Permutation.hs
@@ -33,10 +33,12 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
---import Data.Packed.Matrix hiding(toLists)
-import Data.Packed.Development
+import Numeric.LinearAlgebra.Data hiding(find,size)
+import Numeric.LinearAlgebra.Devel
 
+--import Data.Vector.Storable hiding(Vector,mapM_,reverse)
+--import Data.Vector.Storable hiding(Vector,sum,find,mapM_)
+
 --import Numeric.LinearAlgebra.Linear
 
 --import Control.Monad(when)
@@ -59,6 +61,12 @@
 
 -----------------------------------------------------------------------------
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 data Perm
 type PermHandle = Ptr Perm
 -- | A permutation structure
@@ -226,14 +234,14 @@
 permute :: Permutation -> Vector Double -> Vector Double
 permute (P n p) v = unsafePerformIO $ do
                     r <- createVector n
-                    app2 (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute p' vs vp rs rp) vec v vec r "permute"
+                    (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute p' vs vp rs rp) # v # r #| "permute"
                     return r
 
 -- | apply the inverse permutation to a vector
 inverse_permute :: Permutation -> Vector Double -> Vector Double
 inverse_permute (P n p) v = unsafePerformIO $ do
                     r <- createVector n
-                    app2 (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute_inverse p' vs vp rs rp) vec v vec r "permute"
+                    (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute_inverse p' vs vp rs rp) # v # r #| "permute"
                     return r
 
 -- | multiply two permutations, P = PA * PB
diff --git a/lib/Numeric/GSL/Sort.hs b/lib/Numeric/GSL/Sort.hs
--- a/lib/Numeric/GSL/Sort.hs
+++ b/lib/Numeric/GSL/Sort.hs
@@ -1,7 +1,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Sort
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2015
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -19,19 +19,23 @@
                          ) where
 
 
-import Data.Packed.Vector
-import Data.Packed.Development
+import Numeric.LinearAlgebra.Data
+import Numeric.LinearAlgebra.Devel
 
 import Foreign
 import Foreign.C.Types(CInt(..))
 
 import System.IO.Unsafe(unsafePerformIO)
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
 -- | sort the elements of a vector into ascending order
 sort :: Vector Double -> Vector Double
 sort v = unsafePerformIO $ do
-         r <- createVector (dim v)
-         app2 sort_sort vec v vec r "sort"
+         r <- createVector (size v)
+         sort_sort # v # r #| "sort"
          return r
 
 foreign import ccall "sort-aux.h sort" sort_sort :: CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
diff --git a/lib/Numeric/GSL/Statistics.hs b/lib/Numeric/GSL/Statistics.hs
--- a/lib/Numeric/GSL/Statistics.hs
+++ b/lib/Numeric/GSL/Statistics.hs
@@ -40,10 +40,8 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
---import Data.Packed(Container(..))
-
-import Data.Packed.Development
+import Numeric.LinearAlgebra.Data
+import Numeric.LinearAlgebra.Devel
 
 --import Numeric.GSL.Vector
 --import Numeric.LinearAlgebra.Instances()
@@ -57,19 +55,25 @@
 
 -----------------------------------------------------------------------------
 
+infixl 1 #
+a # b = applyRaw a b
+{-# INLINE (#) #-}
+
+-----------------------------------------------------------------------------
+
 type PD = Ptr Double
 
 -----------------------------------------------------------------------------
 
 getD1 f s v = unsafePerformIO $ do
                 alloca $ \r -> do
-                   app1 (f r) vec v s
+                   (f r) # v #| s
                    r' <- peek r
                    return r'
 
 getD2 f s v w = unsafePerformIO $ do
                    alloca $ \r -> do
-                      app2 (f r) vec v vec w s
+                      (f r) # v # w #| s
                       r' <- peek r
                       return r'