diff --git a/hmatrix-gsl.cabal b/hmatrix-gsl.cabal
--- a/hmatrix-gsl.cabal
+++ b/hmatrix-gsl.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-gsl
-Version:            0.17.0.0
+Version:            0.18.0.1
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -25,12 +25,11 @@
 
 library
 
-    Build-Depends:      base<5, hmatrix>=0.17, array, vector,
+    Build-Depends:      base<5, hmatrix>=0.18, array, vector,
                         process, random
 
 
-    Extensions:         ForeignFunctionInterface,
-                        CPP
+    Extensions:         ForeignFunctionInterface
 
     hs-source-dirs:     src
     Exposed-modules:    Numeric.GSL.Differentiation,
diff --git a/src/Numeric/GSL/Fitting.hs b/src/Numeric/GSL/Fitting.hs
--- a/src/Numeric/GSL/Fitting.hs
+++ b/src/Numeric/GSL/Fitting.hs
@@ -87,7 +87,7 @@
     fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
     jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
     rawpath <- createMatrix RowMajor maxit (2+p)
-    c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n) # xiv # rawpath #|"c_nlfit"
+    (xiv `applyRaw` (rawpath `applyRaw` id)) (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) #|"c_nlfit"
     let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toRows $ dropRows (it-1) path
diff --git a/src/Numeric/GSL/Fourier.hs b/src/Numeric/GSL/Fourier.hs
--- a/src/Numeric/GSL/Fourier.hs
+++ b/src/Numeric/GSL/Fourier.hs
@@ -25,7 +25,7 @@
 
 genfft code v = unsafePerformIO $ do
     r <- createVector (size v)
-    c_fft code # v # r #|"fft"
+    (v `applyRaw` (r `applyRaw` id)) (c_fft code) #|"fft"
     return r
 
 foreign import ccall unsafe "gsl-aux.h fft" c_fft ::  CInt -> TCV (TCV Res)
diff --git a/src/Numeric/GSL/Internal.hs b/src/Numeric/GSL/Internal.hs
--- a/src/Numeric/GSL/Internal.hs
+++ b/src/Numeric/GSL/Internal.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE FlexibleContexts #-}
+
 -- |
 -- Module      :  Numeric.GSL.Internal
 -- Copyright   :  (c) Alberto Ruiz 2009
@@ -23,7 +25,7 @@
     createV,
     createMIO,
     module Numeric.LinearAlgebra.Devel,
-    check,(#),vec, ww2,
+    check,(#),(#!),vec, ww2,
     Res,TV,TM,TCV,TCM
 ) where
 
@@ -86,12 +88,12 @@
 
 createV n fun msg = unsafePerformIO $ do
     r <- createVector n
-    fun # r #| msg
+    (r # id) fun #| msg
     return r
 
 createMIO r c fun msg = do
     res <- createMatrix RowMajor r c
-    fun # res #| msg
+    (res # id) fun #| msg
     return res
 
 --------------------------------------------------------------------------------
@@ -134,4 +136,11 @@
 infixl 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
+
+--infixr 1 #
+--a # b = apply a b
+--{-# INLINE (#) #-}
+
+a #! b = a # b # id
+{-# INLINE (#!) #-}
 
diff --git a/src/Numeric/GSL/Interpolation.hs b/src/Numeric/GSL/Interpolation.hs
--- a/src/Numeric/GSL/Interpolation.hs
+++ b/src/Numeric/GSL/Interpolation.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE MagicHash, UnboxedTuples #-}
+
 {- |
 Module      :  Numeric.GSL.Interpolation
 Copyright   :  (c) Matthew Peddie 2015
@@ -40,6 +42,10 @@
 import Numeric.GSL.Internal
 import System.IO.Unsafe(unsafePerformIO)
 
+-- FIXME
+import qualified Data.Vector.Storable as S
+import GHC.Base (IO(..), realWorld#)
+
 data InterpolationMethod = Linear
                          | Polynomial
                          | CSpline
@@ -58,6 +64,12 @@
 
 dim :: Numeric t => Vector t -> Int
 dim = size
+
+-- FIXME
+appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
+
+unsafeInlinePerformIO (IO f) = case f realWorld# of
+    (# _, x #) -> x
 
 applyCFun hsname cname fun mth xs ys x
   | dim xs /= dim ys = error $
diff --git a/src/Numeric/GSL/LinearAlgebra.hs b/src/Numeric/GSL/LinearAlgebra.hs
--- a/src/Numeric/GSL/LinearAlgebra.hs
+++ b/src/Numeric/GSL/LinearAlgebra.hs
@@ -40,7 +40,7 @@
              -> Vector Double
 randomVector seed dist n = unsafePerformIO $ do
     r <- createVector n
-    c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"
+    (r `applyRaw` id) (c_random_vector (fi seed) ((fi.fromEnum) dist))  #|"randomVector"
     return r
 
 foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
@@ -56,7 +56,7 @@
     charname <- newCString filename
     charfmt <- newCString fmt
     let o = if orderOf m == RowMajor then 1 else 0
-    matrix_fprintf charname charfmt o # m #|"matrix_fprintf"
+    (m `applyRaw` id) (matrix_fprintf charname charfmt o)  #|"matrix_fprintf"
     free charname
     free charfmt
 
@@ -69,7 +69,7 @@
 fscanfVector filename n = do
     charname <- newCString filename
     res <- createVector n
-    gsl_vector_fscanf charname # res #|"gsl_vector_fscanf"
+    (res `applyRaw` id) (gsl_vector_fscanf charname) #|"gsl_vector_fscanf"
     free charname
     return res
 
@@ -80,7 +80,7 @@
 fprintfVector filename fmt v = do
     charname <- newCString filename
     charfmt <- newCString fmt
-    gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf"
+    (v `applyRaw` id) (gsl_vector_fprintf charname charfmt) #|"gsl_vector_fprintf"
     free charname
     free charfmt
 
@@ -91,7 +91,7 @@
 freadVector filename n = do
     charname <- newCString filename
     res <- createVector n
-    gsl_vector_fread charname # res #| "gsl_vector_fread"
+    (res `applyRaw` id) (gsl_vector_fread charname) #| "gsl_vector_fread"
     free charname
     return res
 
@@ -101,7 +101,7 @@
 fwriteVector :: FilePath -> Vector Double -> IO ()
 fwriteVector filename v = do
     charname <- newCString filename
-    gsl_vector_fwrite charname # v #|"gsl_vector_fwrite"
+    (v `applyRaw` id) (gsl_vector_fwrite charname) #|"gsl_vector_fwrite"
     free charname
 
 foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
diff --git a/src/Numeric/GSL/ODE.hs b/src/Numeric/GSL/ODE.hs
--- a/src/Numeric/GSL/ODE.hs
+++ b/src/Numeric/GSL/ODE.hs
@@ -17,7 +17,7 @@
 @
 import Numeric.GSL.ODE
 import Numeric.LinearAlgebra
-import Numeric.LinearAlgebra.Util(mplot)
+import Graphics.Plot(mplot)
 
 xdot t [x,v] = [v, -0.95*x - 0.1*v]
 
diff --git a/src/Numeric/GSL/Polynomials.hs b/src/Numeric/GSL/Polynomials.hs
--- a/src/Numeric/GSL/Polynomials.hs
+++ b/src/Numeric/GSL/Polynomials.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP #-}
 {- |
 Module      :  Numeric.GSL.Polynomials
 Copyright   :  (c) Alberto Ruiz 2006
@@ -24,7 +25,7 @@
 import Foreign.C.Types (CInt(..))
 #endif
 
-{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. 
+{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/.
 
 For example, the three solutions of x^3 + 8 = 0
 
@@ -41,14 +42,14 @@
 0.30901699437494756 :+ (-0.9510565162951535),
 1.0000000000000002 :+ 0.0]
 
--}  
+-}
 polySolve :: [Double] -> [Complex Double]
 polySolve = toList . polySolve' . fromList
 
 polySolve' :: Vector Double -> Vector (Complex Double)
 polySolve' v | size v > 1 = unsafePerformIO $ do
     r <- createVector (size v-1)
-    c_polySolve # v # r #| "polySolve"
+    (v `applyRaw` (r `applyRaw` id)) c_polySolve #| "polySolve"
     return r
              | otherwise = error "polySolve on a polynomial of degree zero"
 
diff --git a/src/Numeric/GSL/Random.hs b/src/Numeric/GSL/Random.hs
--- a/src/Numeric/GSL/Random.hs
+++ b/src/Numeric/GSL/Random.hs
@@ -39,13 +39,13 @@
 gaussianSample :: Seed
                -> Int -- ^ number of rows
                -> Vector Double -- ^ mean vector
-               -> Matrix Double -- ^ covariance matrix
+               -> Herm   Double -- ^ covariance matrix
                -> Matrix Double -- ^ result
 gaussianSample seed n med cov = m where
     c = size med
     meds = konst 1 n `outer` med
     rs = reshape c $ randomVector seed Gaussian (c * n)
-    m = rs <> cholSH cov + meds
+    m = rs <> chol cov + meds
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- uniform distribution.
@@ -87,4 +87,3 @@
 -}
 randn :: Int -> Int -> IO (Matrix Double)
 randn = randm Gaussian
-
diff --git a/src/Numeric/GSL/Vector.hs b/src/Numeric/GSL/Vector.hs
--- a/src/Numeric/GSL/Vector.hs
+++ b/src/Numeric/GSL/Vector.hs
@@ -34,7 +34,7 @@
              -> Vector Double
 randomVector seed dist n = unsafePerformIO $ do
     r <- createVector n
-    c_random_vector_GSL (fi seed) ((fi.fromEnum) dist) # r #|"randomVectorGSL"
+    (r `applyRaw` id) (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) #|"randomVectorGSL"
     return r
 
 foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
@@ -50,7 +50,7 @@
     charname <- newCString filename
     charfmt <- newCString fmt
     let o = if orderOf m == RowMajor then 1 else 0
-    matrix_fprintf charname charfmt o # m #|"matrix_fprintf"
+    (m `applyRaw` id) (matrix_fprintf charname charfmt o) #|"matrix_fprintf"
     free charname
     free charfmt
 
@@ -63,7 +63,7 @@
 fscanfVector filename n = do
     charname <- newCString filename
     res <- createVector n
-    gsl_vector_fscanf charname # res #|"gsl_vector_fscanf"
+    (res `applyRaw` id) (gsl_vector_fscanf charname) #|"gsl_vector_fscanf"
     free charname
     return res
 
@@ -74,7 +74,7 @@
 fprintfVector filename fmt v = do
     charname <- newCString filename
     charfmt <- newCString fmt
-    gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf"
+    (v `applyRaw` id) (gsl_vector_fprintf charname charfmt) #|"gsl_vector_fprintf"
     free charname
     free charfmt
 
@@ -85,7 +85,7 @@
 freadVector filename n = do
     charname <- newCString filename
     res <- createVector n
-    gsl_vector_fread charname # res #|"gsl_vector_fread"
+    (res `applyRaw` id) (gsl_vector_fread charname) #|"gsl_vector_fread"
     free charname
     return res
 
@@ -95,7 +95,7 @@
 fwriteVector :: FilePath -> Vector Double -> IO ()
 fwriteVector filename v = do
     charname <- newCString filename
-    gsl_vector_fwrite charname # v #|"gsl_vector_fwrite"
+    (v `applyRaw` id) (gsl_vector_fwrite charname) #|"gsl_vector_fwrite"
     free charname
 
 foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV