5 files changed, 54 insertions, 75 deletions

diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
index 7ba15aa..6490b25 100644
--- a/examples/benchmarks.hs
+++ b/examples/benchmarks.hs
@@ -18,7 +18,9 @@ time act = do
 
 --------------------------------------------------------------------------------
 
-main = sequence_ [bench1,bench2,bench3,bench4,bench5 1000000 3]
+main = sequence_ [bench1,bench2,bench3,bench4,
+                  bench5 1000000 3,
+                  bench5 100000 50]
 
 w :: Vector Double
 w = constant 1 5000000
@@ -131,19 +133,18 @@ bench4 = do
 --------------------------------------------------------------------------------
 
 op1 a b = a <> trans b
-
-op2 a b = a + trans b
+op2 a b = a  + trans b
 
 timep = time . print . vectorMax . flatten
 
 bench5 n d = do
     putStrLn "-------------------------------------------------------"
-    putStrLn "transpose in multiply"
-    let ms = replicate n ((ident d :: Matrix Double))
-    let mz = replicate n (diag (constant (0::Double) d))
-    timep $ foldl1' (<>) ms
-    timep $ foldl1' op1  ms
-    putStrLn "-------------------------------------------------------"
     putStrLn "transpose in add"
+    let ms = replicate n ((ident d :: Matrix Double))
     timep $ foldl1' (+)  ms
     timep $ foldl1' op2  ms
+    putStrLn "-------------------------------------------------------"
+    putStrLn "transpose in multiply"
+
+    timep $ foldl1' (<>) ms
+    timep $ foldl1' op1  ms
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 3307870..638358b 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -13,7 +13,7 @@ Description:        A purely functional interface to basic linear algebra comput
                     .
                     More information: <http://www.hmatrix.googlepages.com>
 Category:           Math
-tested-with:        GHC ==6.10.0
+tested-with:        GHC ==6.10.1
 
 cabal-version:      >=1.2
 build-type:         Simple
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 477b453..a32a782 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -1,5 +1,5 @@
 {-# OPTIONS_GHC -fglasgow-exts #-}
-{-# LANGUAGE CPP #-}
+{-# LANGUAGE CPP, BangPatterns #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Matrix
@@ -71,6 +71,12 @@ data Matrix t = MC { rows :: {-# UNPACK #-} !Int
 -- MC: preferred by C, fdat may require a transposition
 -- MF: preferred by LAPACK, cdat may require a transposition
 
+xdat MC {cdat = d } = d
+xdat MF {fdat = d } = d
+
+orderOf MF{} = ColumnMajor
+orderOf MC{} = RowMajor
+
 -- | Matrix transpose.
 trans :: Matrix t -> Matrix t
 trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d }
@@ -84,16 +90,10 @@ fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transda
 
 mat = withMatrix
 
-withMatrix MC {rows = r, cols = c, cdat = d } f =
-    withForeignPtr (fptr d) $ \p -> do
-        let m g = do
-            g (fi r) (fi c) p
-        f m
-
-withMatrix MF {rows = r, cols = c, fdat = d } f =
-    withForeignPtr (fptr d) $ \p -> do
+withMatrix a f =
+    withForeignPtr (fptr (xdat a)) $ \p -> do
         let m g = do
-            g (fi r) (fi c) p
+            g (fi (rows a)) (fi (cols a)) p
         f m
 
 {- | Creates a vector by concatenation of rows
@@ -156,6 +156,11 @@ MF {rows = r, cols = c, fdat = v} @@> (i,j)
     | otherwise = v `at` (j*r+i)
 {-# INLINE (@@>) #-}
 
+-- | Unsafe matrix access without range checking
+atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j)
+atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i)
+{-# INLINE atM' #-}
+
 ------------------------------------------------------------------
 
 matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v }
@@ -204,56 +209,52 @@ liftMatrix2 f m1 m2
 compat :: Matrix a -> Matrix b -> Bool
 compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
 
-----------------------------------------------------------------
+------------------------------------------------------------------
 
--- | Optimized matrix computations are provided for elements in the Element class.
+-- | Auxiliary class.
 class (Storable a, Floating a) => Element a where
-    transdata :: Int -> Vector a -> Int -> Vector a
     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
                -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
                -> Matrix a -> Matrix a
+    transdata :: Int -> Vector a -> Int -> Vector a
 
 instance Element Double where
-    transdata = transdataR
     subMatrixD = subMatrixR
+    transdata  = transdata'
 
 instance Element (Complex Double) where
-    transdata  = transdataC
     subMatrixD = subMatrixC
-
-------------------------------------------------------------------
-
-(>|<) :: (Element a) => Int -> Int -> [a] -> Matrix a
-r >|< c = f where
-    f l | dim v == r*c = matrixFromVector ColumnMajor c v
-        | otherwise    = error $ "inconsistent list size = "
-                                 ++show (dim v) ++" in ("++show r++"><"++show c++")"
-        where v = fromList l
+    transdata  = transdata'
 
 -------------------------------------------------------------------
 
-transdataR :: Int -> Vector Double -> Int -> Vector Double
-transdataR = transdataAux ctransR
-
-transdataC :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double)
-transdataC = transdataAux ctransC
-
-transdataAux fun c1 d c2 =
+transdata' :: Storable a => Int -> Vector a -> Int -> Vector a
+transdata' c1 v c2 =
     if noneed
-        then d
+        then v
         else unsafePerformIO $ do
-            v <- createVector (dim d)
-            withForeignPtr (fptr d) $ \pd ->
-                withForeignPtr (fptr v) $ \pv ->
-                    fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
-            -- putStrLn $ "---> transdataAux" ++ show (toList d) ++ show (toList v)
-            return v
-  where r1 = dim d `div` c1
-        r2 = dim d `div` c2
+                w <- createVector (r2*c2)
+                let p = unsafeForeignPtrToPtr (fptr v)
+                    q = unsafeForeignPtrToPtr (fptr w)
+                    go (-1) _ = return ()
+                    go !i (-1) = go (i-1) (c1-1)
+                    go !i !j = do x <- peekElemOff p (i*c1+j)
+                                  pokeElemOff      q (j*c2+i) x
+                                  go i (j-1)
+                go (r1-1) (c1-1)
+                touchForeignPtr (fptr w)
+                return w
+  where r1 = dim v `div` c1
+        r2 = dim v `div` c2
         noneed = r1 == 1 || c1 == 1
 
-foreign import ccall "auxi.h transR" ctransR :: TMM
-foreign import ccall "auxi.h transC" ctransC :: TCMCM
+
+-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int ->  Vector Double #-}
+-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}
+
+-- I don't know how to specialize...
+-- The above pragmas only seem to work on top level defs
+-- Fortunately everything seems to work using the above class
 
 ----------------------------------------------------------------------
 
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index 48b05e8..562c804 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -51,26 +51,6 @@
 #define MEM      2002
 #define BAD_FILE 2003
 
-int transR(KRMAT(x),RMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transR");
-    KDMVIEW(x);
-    DMVIEW(t);
-    int res = gsl_matrix_transpose_memcpy(M(t),M(x));
-    CHECK(res,res);
-    OK
-}
-
-int transC(KCMAT(x),CMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transC");
-    KCMVIEW(x);
-    CMVIEW(t);
-    int res = gsl_matrix_complex_transpose_memcpy(M(t),M(x));
-    CHECK(res,res);
-    OK
-}
-
 
 int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)) {
     REQUIRES(0<=r1 && r1<=r2 && r2<xr && 0<=c1 && c1<=c2 && c2<xc &&
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
index 4698696..e5370b4 100644
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -10,9 +10,6 @@
 #define KCVEC(A) int A##n, const gsl_complex*A##p
 #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
 
-int transR(KRMAT(x),RMAT(t));
-int transC(KCMAT(x),CMAT(t));
-
 int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));
 
 const char * gsl_strerror (const int gsl_errno);