Unpack Matrix type, and add -funsafe flag

-funsafe optionally compiles out the bounds checks on indexing matrices
and vectors. Yields good speedups on tight loops. Not enabled by default.

3 files changed, 43 insertions, 11 deletions

diff --git a/hmatrix.cabal b/hmatrix.cabal
index 2633980..f39b862 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -25,6 +25,10 @@ flag mkl
     description:    Link with Intel's MKL optimized libraries.
     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,7 +40,8 @@ library
         ghc-options:    -fvia-C
 
     Build-Depends:      haskell98
-    Extensions:         ForeignFunctionInterface
+    Extensions:         ForeignFunctionInterface,
+                        CPP
 
     hs-source-dirs:     lib
     Exposed-modules:    Data.Packed,
@@ -107,3 +112,7 @@ library
 
     cc-options:         -O4
     ghc-prof-options:   -auto-all
+
+    if flag(unsafe)
+        cpp-options: -DUNSAFE
+
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index c6129c6..7a35a61 100644
--- 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,8 +62,13 @@ import Data.List(transpose)
 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 @@ withMatrix MF {rows = r, cols = c, fdat = d } f =
 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 @@ fromColumns m = trans . fromRows $ m
 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 @@ infixl 9 @@>
 --    | 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 (@@>) #-}
 
 ------------------------------------------------------------------
 
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 8723367..6274e48 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -100,10 +100,24 @@ n |> l = if length l == n then fromList l else error "|> with wrong size"
 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