diff --git a/README b/README
--- a/README
+++ b/README
@@ -198,8 +198,8 @@
 
 ACKNOWLEDGEMENTS -----------------------------------------------------
 
-I thank Henning Thielemann and all the people in the Haskell mailing
-lists for their help.
+I thank Don Stewart, Henning Thielemann, Bulat Ziganshin and all the people
+in the Haskell mailing lists for their help.
 
 - Nico Mahlo discovered a bug in the eigendecomposition wrapper.
 
@@ -221,3 +221,10 @@
   which are wrong for this architecture.
 
 - Jason Schroeder reported an error in the documentation.
+
+- Bulat Ziganshin gave invaluable help for the ST monad interface to
+  in-place modifications.
+
+- Don Stewart fixed the implementation of the internal data structures
+  to achieve excellent, C-like performance in Haskell functions which
+  explicitly work with the elements of vectors and matrices.
diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
new file mode 100644
--- /dev/null
+++ b/examples/benchmarks.hs
@@ -0,0 +1,103 @@
+{-# OPTIONS -fbang-patterns #-}
+
+-- compile as:
+-- ghc --make -O2 -optc-O2 -fvia-C benchmarks.hs
+-- ghc --make -O benchmarks.hs
+
+import Numeric.LinearAlgebra
+import System.Time
+import System.CPUTime
+import Text.Printf
+import Data.List(foldl1')
+
+
+time act = do
+    t0 <- getCPUTime
+    act
+    t1 <- getCPUTime
+    printf "%.3f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()
+
+--------------------------------------------------------------------------------
+
+main = sequence_ [bench1,bench2,bench3]
+
+w :: Vector Double
+w = constant 1 30000000
+
+bench1 = do
+    putStrLn "Sum of a vector with 30M doubles:"
+    print$ vectorMax w -- evaluate it
+    time $ printf "     BLAS: %.2f: " $ sumVB w
+    time $ printf "  Haskell: %.2f: " $ sumVH w
+    time $ printf "     BLAS: %.2f: " $ sumVB w
+    time $ printf "  Haskell: %.2f: " $ sumVH w
+
+sumVB v = constant 1 (dim v) <.> v
+
+sumVH v = go (d - 1) 0
+     where
+       d = dim v
+       go :: Int -> Double -> Double
+       go 0 s = s + (v @> 0)
+       go !j !s = go (j - 1) (s + (v @> j))
+
+--------------------------------------------------------------------------------
+
+bench2 = do
+    putStrLn "-------------------------------------------------------"
+    putStrLn "Multiplication of 1M different 3x3 matrices:"
+--     putStrLn "from [[]]"
+--     time $ print $ fun (10^6) rot'
+--     putStrLn "from []"
+    time $ print $ manymult (10^6) rot
+    print $ cos (10^6/2)
+
+
+rot' :: Double -> Matrix Double
+rot' a = matrix [[ c,0,s],
+                 [ 0,1,0],
+                 [-s,0,c]]
+    where c = cos a
+          s = sin a
+          matrix = fromLists
+
+rot :: Double -> Matrix Double
+rot a = (3><3) [ c,0,s
+               , 0,1,0
+               ,-s,0,c ]
+    where c = cos a
+          s = sin a
+
+manymult n r = foldl1' (<>) (map r angles)
+    where angles = toList $ linspace n (0,1)
+          -- angles = map (k*) [0..n']
+          -- n' = fromIntegral n - 1
+          -- k  = recip n'
+
+--------------------------------------------------------------------------------
+
+bench3 = do
+    putStrLn "-------------------------------------------------------"
+    putStrLn "foldVector"
+    let v = flatten $ ident 500 :: Vector Double
+    print $ vectorMax v  -- evaluate it
+    let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s
+    putStrLn "indices extraction, dim=0.25M:"
+    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)
+    putStrLn "sum, dim=30M:"
+    --time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)
+    time $ print $ foldVector (\k v s -> v k + s) 0.0 w
+    putStrLn "sum, dim=0.25M:"
+    --time $ print $ foldVector (\k v s -> v k + s) 0.0 v
+    time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)
+
+-- foldVector is slower if it is used in two places. (!?)
+-- this does not happen with foldLoop
+
+foldLoop f s d = go (d - 1) s
+     where
+       go 0 s = f (0::Int) s
+       go !j !s = go (j - 1) (f j s)
+
+foldVector f s v = foldLoop g s (dim v)
+    where g !k !s = f k (v@>) s
diff --git a/examples/experiments/parallel.hs b/examples/experiments/parallel.hs
--- a/examples/experiments/parallel.hs
+++ b/examples/experiments/parallel.hs
@@ -5,24 +5,32 @@
 
 inParallel = parMap rwhnf id
 
-parMul x y = fromBlocks [inParallel[x <> y1, x <> y2]]
-    where p = cols y `div` 2
-          (y1,y2) = splitColumnsAt p y
+parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]
+    where ys = splitColumns p y
 
+
 main = do
     n <- (read . head) `fmap` getArgs
     let m = ident n :: Matrix Double
-    time $ print $ vectorMax $ takeDiag $ parMul m m
     time $ print $ vectorMax $ takeDiag $ m <> m
-
-a = (2><3) [1..6::Double]
-b = (3><4) [1..12::Double]
-
-splitRowsAt p m    = (takeRows p m, dropRows p m)
-splitColumnsAt p m = (takeColumns p m, dropColumns p m)
+    time $ print $ vectorMax $ takeDiag $ parMul 2 m m
+    time $ print $ vectorMax $ takeDiag $ parMul 4 m m
+    time $ print $ vectorMax $ takeDiag $ parMul 8 m m
 
 time act = do
     t0 <- getClockTime
     act
     t1 <- getClockTime
     print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0
+
+splitColumns n m = splitColumns' (f n (cols m)) m
+    where
+    splitColumns' [] m = []
+    splitColumns' ((a,b):rest) m = subMatrix (0,a) (rows m, b-a+1) m : splitColumns' rest m
+
+    f :: Int -> Int -> [(Int,Int)]
+    f n c = zip ks (map pred $ tail ks)
+        where ks = map round $ toList $ linspace (fromIntegral n+1) (0,fromIntegral c)
+
+splitRowsAt p m    = (takeRows p m, dropRows p m)
+splitColumnsAt p m = (takeColumns p m, dropColumns p m)
diff --git a/examples/inplace.hs b/examples/inplace.hs
new file mode 100644
--- /dev/null
+++ b/examples/inplace.hs
@@ -0,0 +1,152 @@
+-- some tests of the interface for pure
+-- computations with inplace updates
+
+import Numeric.LinearAlgebra
+import Data.Packed.ST
+import Data.Packed.Convert
+
+import Data.Array.Unboxed
+import Data.Array.ST
+import Control.Monad.ST
+import Control.Monad
+
+main = sequence_[
+    print test1,
+    print test2,
+    print test3,
+    print test4,
+    test5,
+    test6,
+    print test7,
+    test8,
+    test0]
+
+-- helper functions
+vector l = fromList l :: Vector Double
+norm v = pnorm PNorm2 v
+
+-- hmatrix vector and matrix
+v = vector [1..10]
+m = (5><10) [1..50::Double]
+
+----------------------------------------------------------------------
+
+-- vector creation by in-place updates on a copy of the argument
+test1 = fun v
+
+fun :: Element t => Vector t -> Vector t
+fun x = runSTVector $ do
+    a <- thawVector x
+    mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]
+    return a
+
+-- another example: creation of an antidiagonal matrix from a list
+test2 = antiDiag 5 8 [1..] :: Matrix Double
+
+antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
+antiDiag r c l = runSTMatrix $ do
+    m <- newMatrix 0 r c
+    let d = min r c - 1
+    sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l
+    return m
+
+-- using vector or matrix functions on mutable objects requires freezing:
+test3 = g1 v
+
+g1 x = runST $ do
+    a <- thawVector x
+    writeVector a (dim x -1) 0
+    b <- freezeVector a
+    return (norm b)
+
+--  another possibility:
+test4 = g2 v
+
+g2 x = runST $ do
+    a <- thawVector x
+    writeVector a (dim x -1) 0
+    t <- liftSTVector norm a
+    return t
+
+--------------------------------------------------------------
+
+-- haskell arrays
+hv = listArray (0,9) [1..10::Double]
+hm = listArray ((0,0),(4,9)) [1..50::Double]
+
+
+
+-- conversion from standard Haskell arrays
+test5 = do
+    print $ norm (vectorFromArray hv)
+    print $ norm v
+    print $ rcond (matrixFromArray hm)
+    print $ rcond m
+
+
+-- conversion to mutable ST arrays
+test6 = do
+    let y = clearColumn m 1
+    print y
+    print (matrixFromArray y)
+
+clearColumn x c = runSTUArray $ do
+    a <- mArrayFromMatrix x
+    forM_ [0..rows x-1] $ \i->
+        writeArray a (i,c) (0::Double)
+    return a
+
+-- hmatrix functions applied to mutable ST arrays
+test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)
+
+unitary v = runSTUArray $ do
+    a <- thaw v
+    n <- norm `fmap` vectorFromMArray a
+    b <- mapArray (/n) a
+    return b
+
+-------------------------------------------------
+
+-- (just to check that they are not affected)
+test0 = do
+    print v
+    print m
+    --print hv
+    --print hm
+
+-------------------------------------------------
+
+histogram n ds = runSTVector $ do
+    h <- newVector (0::Double) n -- number of bins
+    let inc k = modifyVector h k (+1)
+    mapM_ inc ds
+    return h
+
+-- check that newVector is really called with a fresh new array
+histoCheck ds = runSTVector $ do
+    h <- newVector (0::Double) 15 -- > constant for this test
+    let inc k = modifyVector h k (+1)
+    mapM_ inc ds
+    return h
+
+hc = fromList [1 .. 15::Double]
+
+-- check that thawVector creates a new array
+histoCheck2 ds = runSTVector $ do
+    h <- thawVector hc
+    let inc k = modifyVector h k (+1)
+    mapM_ inc ds
+    return h
+
+test8 = do
+    let ds = [0..14]
+    print $ histogram 15 ds
+    print $ histogram 15 ds
+    print $ histogram 15 ds
+    print $ histoCheck ds
+    print $ histoCheck ds
+    print $ histoCheck ds
+    print $ histoCheck2 ds
+    print $ histoCheck2 ds
+    print $ histoCheck2 ds
+    putStrLn "----------------------"
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,19 +1,19 @@
 Name:               hmatrix
-Version:            0.3.0.0
+Version:            0.4.0.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
 Maintainer:         Alberto Ruiz <aruiz@um.es>
 Stability:          provisional
-Homepage:           http://alberrto.googlepages.com/gslhaskell
+Homepage:           http://www.hmatrix.googlepages.com
 Synopsis:           Linear algebra and numerical computations
 Description:        A purely functional interface to basic linear algebra computations
                     and other numerical routines, internally implemented using
                     GSL, BLAS and LAPACK.
                     .
-                    More information: <http://alberrto.googlepages.com/gslhaskell>
+                    More information: <http://www.hmatrix.googlepages.com>
 Category:           Numerical, Math
-tested-with:        GHC ==6.8.2
+tested-with:        GHC ==6.8.3
 
 cabal-version:      >=1.2
 build-type:         Simple
@@ -25,6 +25,14 @@
     description:    Link with Intel's MKL optimized libraries.
     default:        False
 
+flag gsl
+    description:    Link with GSL unoptimized blas.
+    default:        False
+
+flag unsafe
+    description: Compile the library with bound checking disabled.
+    default: False
+
 library
     if flag(splitBase)
       build-depends:    base >= 3, array, QuickCheck, HUnit, storable-complex
@@ -36,19 +44,19 @@
         ghc-options:    -fvia-C
 
     Build-Depends:      haskell98
-    Extensions:         ForeignFunctionInterface
+    Extensions:         ForeignFunctionInterface,
+                        CPP
 
     hs-source-dirs:     lib
     Exposed-modules:    Data.Packed,
                         Data.Packed.Vector,
                         Data.Packed.Matrix,
-                        Numeric.GSL.Vector,
-                        Numeric.GSL.Matrix,
                         Numeric.GSL.Differentiation,
                         Numeric.GSL.Integration,
                         Numeric.GSL.Fourier,
                         Numeric.GSL.Polynomials,
                         Numeric.GSL.Minimization,
+                        Numeric.GSL.Vector,
                         Numeric.GSL.Special,
                         Numeric.GSL.Special.Gamma,
                         Numeric.GSL.Special.Erf,
@@ -85,12 +93,15 @@
                         Numeric.LinearAlgebra.Algorithms,
                         Graphics.Plot,
                         Numeric.LinearAlgebra.Tests
+                        Data.Packed.Convert
+                        Data.Packed.ST
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Vector,
                         Data.Packed.Internal.Matrix,
                         Numeric.GSL.Special.Internal,
-                        Numeric.LinearAlgebra.Tests.Instances
+                        Numeric.GSL.Matrix,
+                        Numeric.LinearAlgebra.Tests.Instances,
                         Numeric.LinearAlgebra.Tests.Properties
     C-sources:          lib/Data/Packed/Internal/auxi.c,
                         lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
@@ -101,7 +112,14 @@
       else
         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core
     else
-      extra-libraries:  gsl blas lapack
+      if flag(gsl)
+        extra-libraries:  gsl gslcblas lapack
+      else
+        extra-libraries:  gsl blas lapack
 
     cc-options:         -O4
     ghc-prof-options:   -auto-all
+
+    if flag(unsafe)
+        cpp-options: -DUNSAFE
+
diff --git a/lib/Data/Packed/Convert.hs b/lib/Data/Packed/Convert.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Packed/Convert.hs
@@ -0,0 +1,102 @@
+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Convert
+-- Copyright   :  (c) Alberto Ruiz 2008
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Conversion of Vectors and Matrices to and from the standard Haskell arrays.
+-- (provisional)
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Convert (
+    arrayFromVector, vectorFromArray,
+    mArrayFromVector, vectorFromMArray,
+    vectorToStorableArray, storableArrayToVector,
+    arrayFromMatrix, matrixFromArray,
+    mArrayFromMatrix, matrixFromMArray,
+--    matrixToStorableArray, storableArrayToMatrix
+) where
+
+import Data.Packed.Internal
+import Data.Array.Storable
+import Foreign
+import Control.Monad.ST
+import Data.Array.ST
+import Data.Array.IArray
+import Data.Array.Unboxed
+
+-- | Creates a StorableArray indexed from 0 to dim -1.
+--   (Memory is efficiently copied, so you can then freely modify the obtained array)
+vectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t)
+vectorToStorableArray v = do
+    r <- cloneVector v
+    unsafeForeignPtrToStorableArray (fptr r) (0,dim r -1)
+
+-- | Creates a Vector from a StorableArray.
+--   (Memory is efficiently copied, so posterior changes in the array will not affect the result)
+storableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t)
+storableArrayToVector s = do
+    (a,b) <- getBounds s
+    let n = (b-a+1)
+    r <- createVector n
+    withStorableArray s $ \p -> do
+        let f _ d =  copyArray d p n >> return 0
+        app1 f vec r "storableArrayToVector"
+    return r
+
+
+unsafeVectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t)
+unsafeVectorToStorableArray v = unsafeForeignPtrToStorableArray (fptr v) (0,dim v -1)
+
+--unsafeStorableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t)
+--unsafeStorableArrayToVector s = undefined -- the foreign ptr of Storable Array is not available?
+
+-----------------------------------------------------------------
+-- provisional, we need Unboxed arrays for Complex Double
+
+
+unsafeFreeze' :: (MArray a e m, Ix i) => a i e -> m (Array i e)
+unsafeFreeze' = unsafeFreeze
+
+-- | creates an immutable Array from an hmatrix Vector (to do: unboxed)
+arrayFromVector :: (Storable t) => Vector t -> Array Int t
+arrayFromVector x = runSTArray (mArrayFromVector x)
+
+-- | creates a mutable array from an hmatrix Vector (to do: unboxed)
+mArrayFromVector :: (MArray b t (ST s), Storable t) => Vector t -> ST s (b Int t)
+mArrayFromVector v = unsafeThaw =<< unsafeIOToST ( unsafeFreeze' =<< (vectorToStorableArray $ v))
+
+
+-- (creates an hmatrix Vector from an immutable array (to do: unboxed))
+vectorFromArray :: (Storable t) => Array Int t -> Vector t
+vectorFromArray a = unsafePerformIO $ storableArrayToVector =<< unsafeThaw a
+
+-- | creates a mutable Array from an hmatrix Vector for manipulation with runSTUArray (to do: unboxed)
+vectorFromMArray :: (Storable t, MArray a t (ST s)) => a Int t -> ST s (Vector t)
+vectorFromMArray x = fmap vectorFromArray (unsafeFreeze' x)
+
+--------------------------------------------------------------------
+-- provisional
+
+matrixFromArray :: UArray (Int, Int) Double -> Matrix Double
+matrixFromArray m = reshape c . fromList . elems $ m
+    where ((r1,c1),(r2,c2)) = bounds m
+          r = r2-r1+1
+          c = c2-c1+1
+
+arrayFromMatrix :: Matrix Double -> UArray (Int, Int) Double
+arrayFromMatrix m = listArray ((0,0),(rows m -1, cols m -1)) (toList $ flatten m)
+
+
+mArrayFromMatrix :: (MArray b Double m) => Matrix Double -> m (b (Int, Int) Double)
+mArrayFromMatrix = unsafeThaw . arrayFromMatrix
+
+matrixFromMArray :: (MArray a Double (ST s)) => a (Int,Int) Double -> ST s (Matrix Double)
+matrixFromMArray x = fmap matrixFromArray (unsafeFreeze x)
diff --git a/lib/Data/Packed/Internal.hs b/lib/Data/Packed/Internal.hs
--- a/lib/Data/Packed/Internal.hs
+++ b/lib/Data/Packed/Internal.hs
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -1,4 +1,5 @@
 {-# OPTIONS_GHC -fglasgow-exts #-}
+{-# LANGUAGE CPP #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Matrix
@@ -61,9 +62,14 @@
 data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
 
 -- | Matrix representation suitable for GSL and LAPACK computations.
-data Matrix t = MC { rows :: Int, cols :: Int, cdat :: Vector t }
-              | MF { rows :: Int, cols :: Int, fdat :: Vector t }
+data Matrix t = MC { rows :: {-# UNPACK #-} !Int
+                   , cols :: {-# UNPACK #-} !Int
+                   , cdat :: {-# UNPACK #-} !(Vector t) }
 
+              | MF { rows :: {-# UNPACK #-} !Int
+                   , cols :: {-# UNPACK #-} !Int
+                   , fdat :: {-# UNPACK #-} !(Vector t) }
+
 -- MC: preferred by C, fdat may require a transposition
 -- MF: preferred by LAPACK, cdat may require a transposition
 
@@ -100,7 +106,6 @@
 flatten :: Element t => Matrix t -> Vector t
 flatten = cdat . cmat
 
-
 type Mt t s = Int -> Int -> Ptr t -> s
 -- not yet admitted by my haddock version
 -- infixr 6 ::>
@@ -133,7 +138,6 @@
 toColumns :: Element t => Matrix t -> [Vector t]
 toColumns m = toRows . trans $ m
 
-
 -- | Reads a matrix position.
 (@@>) :: Storable t => Matrix t -> (Int,Int) -> t
 infixl 9 @@>
@@ -142,12 +146,17 @@
 --    | otherwise   = cdat m `at` (i*c+j)
 
 MC {rows = r, cols = c, cdat = v} @@> (i,j)
-    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
-    | otherwise                  = v `at` (i*c+j)
+    | safe      = if i<0 || i>=r || j<0 || j>=c
+                    then error "matrix indexing out of range"
+                    else v `at` (i*c+j)
+    | otherwise = v `at` (i*c+j)
 
 MF {rows = r, cols = c, fdat = v} @@> (i,j)
-    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
-    | otherwise                  = v `at` (j*r+i)
+    | safe      = if i<0 || i>=r || j<0 || j>=c
+                    then error "matrix indexing out of range"
+                    else v `at` (j*r+i)
+    | otherwise = v `at` (j*r+i)
+{-# INLINE (@@>) #-}
 
 ------------------------------------------------------------------
 
@@ -254,8 +263,8 @@
         r2 = dim d `div` c2
         noneed = r1 == 1 || c1 == 1
 
-foreign import ccall unsafe "auxi.h transR" ctransR :: TMM
-foreign import ccall unsafe "auxi.h transC" ctransC :: TCMCM
+foreign import ccall "auxi.h transR" ctransR :: TMM
+foreign import ccall "auxi.h transC" ctransC :: TCMCM
 
 ------------------------------------------------------------------
 
@@ -275,10 +284,10 @@
     return r
 
 multiplyR = multiplyAux cmultiplyR
-foreign import ccall unsafe "auxi.h multiplyR" cmultiplyR :: TauxMul Double
+foreign import ccall "auxi.h multiplyR" cmultiplyR :: TauxMul Double
 
 multiplyC = multiplyAux cmultiplyC
-foreign import ccall unsafe "auxi.h multiplyC" cmultiplyC :: TauxMul (Complex Double)
+foreign import ccall "auxi.h multiplyC" cmultiplyC :: TauxMul (Complex Double)
 
 -- | matrix product
 multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+{-# LANGUAGE MagicHash, CPP, UnboxedTuples #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Vector
@@ -21,10 +21,20 @@
 import Complex
 import Control.Monad(when)
 
+#if __GLASGOW_HASKELL__ >= 605
+import GHC.ForeignPtr           (mallocPlainForeignPtrBytes)
+#else
+import Foreign.ForeignPtr       (mallocForeignPtrBytes)
+#endif
+
+import GHC.Base
+import GHC.IOBase
+
 -- | A one-dimensional array of objects stored in a contiguous memory block.
-data Vector t = V { dim  :: Int              -- ^ number of elements
-                  , fptr :: ForeignPtr t     -- ^ foreign pointer to the memory block
-                  }
+data Vector t =
+    V { dim  :: {-# UNPACK #-} !Int               -- ^ number of elements
+      , fptr :: {-# UNPACK #-} !(ForeignPtr t)    -- ^ foreign pointer to the memory block
+      }
 
 vec = withVector
 
@@ -37,8 +47,20 @@
 createVector :: Storable a => Int -> IO (Vector a)
 createVector n = do
     when (n <= 0) $ error ("trying to createVector of dim "++show n)
-    fp <- mallocForeignPtrArray n
+    fp <- doMalloc undefined
     return $ V n fp
+  where
+    --
+    -- Use the much cheaper Haskell heap allocated storage
+    -- for foreign pointer space we control
+    --
+    doMalloc :: Storable b => b -> IO (ForeignPtr b)
+    doMalloc dummy = do
+#if __GLASGOW_HASKELL__ >= 605
+        mallocPlainForeignPtrBytes (n * sizeOf dummy)
+#else
+        mallocForeignPtrBytes      (n * sizeOf dummy)
+#endif
 
 {- | creates a Vector from a list:
 
@@ -53,8 +75,13 @@
     app1 f vec v "fromList"
     return v
 
-safeRead v = unsafePerformIO . withForeignPtr (fptr v)
+safeRead v = inlinePerformIO . withForeignPtr (fptr v)
+{-# INLINE safeRead #-}
 
+inlinePerformIO :: IO a -> a
+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
+{-# INLINE inlinePerformIO #-}
+
 {- | extracts the Vector elements to a list
 
 @> toList (linspace 5 (1,10))
@@ -73,10 +100,25 @@
 at' :: Storable a => Vector a -> Int -> a
 at' v n = safeRead v $ flip peekElemOff n
 
+--
+-- turn off bounds checking with -funsafe at configure time.
+-- ghc will optimise away the salways true case at compile time.
+--
+#if defined(UNSAFE)
+safe :: Bool
+safe = False
+#else
+safe = True
+#endif
+
 -- | access to Vector elements with range checking.
 at :: Storable a => Vector a -> Int -> a
-at v n | n >= 0 && n < dim v = at' v n
-       | otherwise          = error "vector index out of range"
+at v n
+    | safe      = if n >= 0 && n < dim v
+                    then at' v n
+                    else error "vector index out of range"
+    | otherwise = at' v n
+{-# INLINE at #-}
 
 {- | takes a number of consecutive elements from a Vector
 
@@ -146,3 +188,10 @@
 liftVector2 f u v = fromList $ zipWith f (toList u) (toList v)
 
 -----------------------------------------------------------------
+
+cloneVector :: Storable t => Vector t -> IO (Vector t)
+cloneVector (v@V {dim=n}) = do
+        r <- createVector n
+        let f _ s _ d =  copyArray d s n >> return 0
+        app2 f vec v vec r "cloneVector"
+        return r
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Packed/ST.hs
@@ -0,0 +1,160 @@
+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.ST
+-- Copyright   :  (c) Alberto Ruiz 2008
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- In-place manipulation inside the ST monad.
+-- See examples/inplace.hs in the distribution.
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.ST (
+    -- * Mutable Vectors
+    STVector, newVector, thawVector, freezeVector, runSTVector,
+    readVector, writeVector, modifyVector, liftSTVector,
+    -- * Mutable Matrices
+    STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
+    readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
+    -- * Unsafe functions
+    unsafeReadVector, unsafeWriteVector,
+    unsafeThawVector, unsafeFreezeVector,
+    unsafeReadMatrix, unsafeWriteMatrix,
+    unsafeThawMatrix, unsafeFreezeMatrix
+) where
+
+import Data.Packed.Internal
+import Data.Array.Storable
+import Control.Monad.ST
+import Data.Array.ST
+import Foreign
+
+{-# INLINE ioReadV #-}
+ioReadV :: Storable t => Vector t -> Int -> IO t
+ioReadV v k = withForeignPtr (fptr v) $ \s -> peekElemOff s k
+
+{-# INLINE ioWriteV #-}
+ioWriteV :: Storable t => Vector t -> Int -> t -> IO ()
+ioWriteV v k x = withForeignPtr (fptr v) $ \s -> pokeElemOff s k x
+
+newtype STVector s t = STVector (Vector t)
+
+thawVector :: Storable t => Vector t -> ST s (STVector s t)
+thawVector = unsafeIOToST . fmap STVector . cloneVector
+
+unsafeThawVector :: Storable t => Vector t -> ST s (STVector s t)
+unsafeThawVector = unsafeIOToST . return . STVector
+
+runSTVector :: Storable t => (forall s . ST s (STVector s t)) -> Vector t
+runSTVector st = runST (st >>= unsafeFreezeVector)
+
+{-# INLINE unsafeReadVector #-}
+unsafeReadVector :: Storable t => STVector s t -> Int -> ST s t
+unsafeReadVector   (STVector x) = unsafeIOToST . ioReadV x
+
+{-# INLINE unsafeWriteVector #-}
+unsafeWriteVector :: Storable t => STVector s t -> Int -> t -> ST s ()
+unsafeWriteVector  (STVector x) k = unsafeIOToST . ioWriteV x k
+
+{-# INLINE modifyVector #-}
+modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s ()
+modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k
+
+liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a
+liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x
+
+freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
+freezeVector v = liftSTVector id v
+
+unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
+unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
+
+{-# INLINE safeIndexV #-}
+safeIndexV f (STVector v) k
+    | k < 0 || k>= dim v = error $ "out of range error in vector (dim="
+                                   ++show (dim v)++", pos="++show k++")"
+    | otherwise = f (STVector v) k
+
+{-# INLINE readVector #-}
+readVector :: Storable t => STVector s t -> Int -> ST s t
+readVector = safeIndexV unsafeReadVector
+
+{-# INLINE writeVector #-}
+writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
+writeVector = safeIndexV unsafeWriteVector
+
+{-# NOINLINE newVector #-}
+newVector :: Element t => t -> Int -> ST s (STVector s t)
+newVector v = unsafeThawVector . constant v
+
+-------------------------------------------------------------------------
+
+{-# INLINE ioReadM #-}
+ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t
+ioReadM (MC nr nc cv) r c = ioReadV cv (r*nc+c)
+ioReadM (MF nr nc fv) r c = ioReadV fv (c*nr+r)
+
+{-# INLINE ioWriteM #-}
+ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO ()
+ioWriteM (MC nr nc cv) r c val = ioWriteV cv (r*nc+c) val
+ioWriteM (MF nr nc fv) r c val = ioWriteV fv (c*nr+r) val
+
+newtype STMatrix s t = STMatrix (Matrix t)
+
+thawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
+thawMatrix = unsafeIOToST . fmap STMatrix . cloneMatrix
+
+unsafeThawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
+unsafeThawMatrix = unsafeIOToST . return . STMatrix
+
+runSTMatrix :: Storable t => (forall s . ST s (STMatrix s t)) -> Matrix t
+runSTMatrix st = runST (st >>= unsafeFreezeMatrix)
+
+{-# INLINE unsafeReadMatrix #-}
+unsafeReadMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
+unsafeReadMatrix   (STMatrix x) r = unsafeIOToST . ioReadM x r
+
+{-# INLINE unsafeWriteMatrix #-}
+unsafeWriteMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
+unsafeWriteMatrix  (STMatrix x) r c = unsafeIOToST . ioWriteM x r c
+
+{-# INLINE modifyMatrix #-}
+modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
+modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c
+
+liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a
+liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x
+
+unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
+unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
+
+freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
+freezeMatrix m = liftSTMatrix id m
+
+cloneMatrix (MC r c d) = cloneVector d >>= return . MC r c
+cloneMatrix (MF r c d) = cloneVector d >>= return . MF r c
+
+{-# INLINE safeIndexM #-}
+safeIndexM f (STMatrix m) r c
+    | r<0 || r>=rows m ||
+      c<0 || c>=cols m = error $ "out of range error in matrix (size="
+                                 ++show (rows m,cols m)++", pos="++show (r,c)++")"
+    | otherwise = f (STMatrix m) r c
+
+{-# INLINE readMatrix #-}
+readMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
+readMatrix = safeIndexM unsafeReadMatrix
+
+{-# INLINE writeMatrix #-}
+writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
+writeMatrix = safeIndexM unsafeWriteMatrix
+
+{-# NOINLINE newMatrix #-}
+newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t)
+newMatrix v r c = unsafeThawMatrix . reshape c . constant v $ r*c
diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs
--- a/lib/Numeric/GSL/Fourier.hs
+++ b/lib/Numeric/GSL/Fourier.hs
@@ -35,7 +35,7 @@
 
 {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
 
-@> fft ('GSL.Matrix.fromList' [1,2,3,4])
+@> fft ('fromList' [1,2,3,4])
 vector (4) [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]@
 
 -}
diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
--- a/lib/Numeric/GSL/Matrix.hs
+++ b/lib/Numeric/GSL/Matrix.hs
@@ -8,7 +8,7 @@
 -- Stability   :  provisional
 -- Portability :  portable (uses FFI)
 --
--- A few linear algebra computations based on the Numeric.GSL (<http://www.gnu.org/software/Numeric.GSL>).
+-- A few linear algebra computations based on GSL.
 --
 -----------------------------------------------------------------------------
 -- #hide
@@ -92,23 +92,6 @@
 
 {- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/:
 
-@\> let (u,s,v) = svdg $ 'fromLists' [[1,2,3],[-4,1,7]]
-\
-\> u
-0.310 -0.951
-0.951  0.310
-\
-\> s
-8.497 2.792
-\
-\> v
--0.411 -0.785
- 0.185 -0.570
- 0.893 -0.243
-\
-\> u \<\> 'diag' s \<\> 'trans' v
- 1. 2. 3.
--4. 1. 7.@
 
 -}
 svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
@@ -128,20 +111,6 @@
 
 {- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/.
 
-@\> let (q,r) = qr $ 'fromLists' [[1,3,5,7],[2,0,-2,4]]
-\
-\> q
--0.447 -0.894
--0.894  0.447
-\
-\> r
--2.236 -1.342 -0.447 -6.708
-    0. -2.683 -5.367 -4.472
-\
-\> q \<\> r
-1.000 3.000  5.000 7.000
-2.000    0. -2.000 4.000@
-
 -}
 qr :: Matrix Double -> (Matrix Double, Matrix Double)
 qr = qr' . cmat
@@ -299,7 +268,7 @@
   2.+3.i   -7.   0.
       1.    2.  -3.
       1.  -1.i  2.i
-\ 
+\ -- CPP
 \> l \<\> u
  2.+3.i   -7.   0.
      1.    2.  -3.
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -31,14 +31,14 @@
 {- | The method of Nelder and Mead, implemented by /gsl_multimin_fminimizer_nmsimplex/. The gradient of the function is not required. This is the example in the GSL manual:
 
 @minimize f xi = minimizeNMSimplex f xi (replicate (length xi) 1) 1e-2 100
-\ 
+\  -- 
 f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
-\ 
+\  --
 main = do
     let (s,p) = minimize f [5,7]
     print s
     print p
-\ 
+\  --
 \> main
 [0.9920430849306285,1.9969168063253164]
 0. 512.500    1.082 6.500    5.
@@ -104,14 +104,14 @@
 
 @minimize = minimizeConjugateGradient 1E-2 1E-4 1E-3 30
 f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
-\ 
+\  --
 df [x,y] = [20*(x-1), 40*(y-2)]
-\  
+\   --
 main = do
     let (s,p) = minimize f df [5,7]
     print s
     print p
-\ 
+\  --
 \> main
 [1.0,2.0]
  0. 687.848 4.996 6.991
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -300,12 +300,12 @@
 @\> let m = 'fromLists' [[1,0,    0]
                     ,[0,1,    0]
                     ,[0,0,1e-10]]
-\ 
+\  --
 \> 'pinv' m 
 1. 0.           0.
 0. 1.           0.
 0. 0. 10000000000.
-\ 
+\  --
 \> pinvTol 1E8 m
 1. 0. 0.
 0. 1. 0.
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -9,7 +9,7 @@
 Stability   :  provisional
 Portability :  portable
 
-This module exports Show, Eq, Num, Fractional, and Floating instances for Vector and Matrix.
+This module exports Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix.
 
 In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.