data structures simplification

3 files changed, 95 insertions, 46 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index fbab33c..010c40b 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -62,21 +62,35 @@ 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, fdat :: Vector t }
-              | MF { rows :: Int, cols :: Int, fdat :: Vector t, cdat :: Vector t }
+data Matrix t = MC { rows :: Int, cols :: Int, cdat :: Vector t }
+              | MF { rows :: Int, cols :: Int, fdat :: Vector t }
 
 -- MC: preferred by C, fdat may require a transposition
 -- MF: preferred by LAPACK, cdat may require a transposition
 
 -- | Matrix transpose.
 trans :: Matrix t -> Matrix t
-trans MC {rows = r, cols = c, cdat = d, fdat = dt } = MF {rows = c, cols = r, fdat = d, cdat = dt }
-trans MF {rows = r, cols = c, fdat = d, cdat = dt } = MC {rows = c, cols = r, cdat = d, fdat = dt }
+trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d }
+trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d }
 
-dat MC { cdat = d } = d
-dat MF { fdat = d } = d
+cmat m@MC{} = m
+cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c}
+
+fmat m@MF{} = m
+fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r}
+
+matc m f = f (rows m) (cols m) (ptr (cdat m))
+matf m f = f (rows m) (cols m) (ptr (fdat m))
+
+
+{- | Creates a vector by concatenation of rows
+
+@\> flatten ('ident' 3)
+9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
+-}
+flatten :: Element t => Matrix t -> Vector t
+flatten = cdat . cmat
 
-mat d m f = f (rows m) (cols m) (ptr (d m))
 
 type Mt t s = Int -> Int -> Ptr t -> s
 -- not yet admitted by my haddock version
@@ -85,7 +99,7 @@ type Mt t s = Int -> Int -> Ptr t -> s
 
 -- | the inverse of 'Data.Packed.Matrix.fromLists'
 toLists :: (Element t) => Matrix t -> [[t]]
-toLists m = partit (cols m) . toList . cdat $ m
+toLists m = partit (cols m) . toList . flatten $ m
 
 -- | creates a Matrix from a list of vectors
 fromRows :: Element t => [Vector t] -> Matrix t
@@ -96,7 +110,7 @@ fromRows vs = case common dim vs of
 -- | extracts the rows of a matrix as a list of vectors
 toRows :: Element t => Matrix t -> [Vector t]
 toRows m = toRows' 0 where
-    v = cdat m
+    v = flatten $ m
     r = rows m
     c = cols m
     toRows' k | k == r*c  = []
@@ -128,12 +142,12 @@ MF {rows = r, cols = c, fdat = v} @@> (i,j)
 
 ------------------------------------------------------------------
 
-matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v, fdat = transdata c v r }
+matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v }
     where (d,m) = dim v `divMod` c
           r | m==0 = d
             | otherwise = error "matrixFromVector"
 
-matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v, cdat = transdata r v c }
+matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v }
     where (d,m) = dim v `divMod` c
           r | m==0 = d
             | otherwise = error "matrixFromVector"
@@ -167,8 +181,8 @@ liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vec
 liftMatrix2 f m1 m2
     | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
     | otherwise = case m1 of
-        MC {} -> matrixFromVector RowMajor    (cols m1) (f (cdat m1) (cdat m2))
-        MF {} -> matrixFromVector ColumnMajor (cols m1) (f (fdat m1) (fdat m2))
+        MC {} -> matrixFromVector RowMajor    (cols m1) (f (cdat m1) (flatten m2))
+        MF {} -> matrixFromVector ColumnMajor (cols m1) (f (fdat m1) ((fdat.fmat) m2))
 
 
 compat :: Matrix a -> Matrix b -> Bool
@@ -222,8 +236,8 @@ transdataAux fun c1 d c2 =
         then d
         else unsafePerformIO $ do
             v <- createVector (dim d)
-            fun r1 c1 (ptr d) r2 c2 (ptr v) // check "transdataAux" [d]
-            --putStrLn "---> transdataAux"
+            fun r1 c1 (ptr d) r2 c2 (ptr v) // check "transdataAux" [d,v]
+            -- putStrLn $ "---> transdataAux" ++ show (toList d) ++ show (toList v)
             return v
   where r1 = dim d `div` c1
         r2 = dim d `div` c2
@@ -239,11 +253,14 @@ foreign import ccall safe "auxi.h transC"
 gmatC MF {rows = r, cols = c, fdat = d} f = f 1 c r (ptr d)
 gmatC MC {rows = r, cols = c, cdat = d} f = f 0 r c (ptr d)
 
+dtt MC { cdat = d } = d
+dtt MF { fdat = d } = d
+
 multiplyAux fun a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++
                                       show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
     r <- createMatrix RowMajor (rows a) (cols b)
-    fun // gmatC a // gmatC b // mat dat r // check "multiplyAux" [dat a, dat b]
+    fun // gmatC a // gmatC b // matc r // check "multiplyAux" [dtt a, dtt b, cdat r]
     return r
 
 multiplyR = multiplyAux cmultiplyR
@@ -273,18 +290,19 @@ multiply = multiplyD
 
 -- | extraction of a submatrix from a real matrix
 subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double
-subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do
+subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do
     r <- createMatrix RowMajor rt ct
-    c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r]
+    let x = cmat x'
+    c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // matc x // matc r // check "subMatrixR" [cdat x]
     return r
 foreign import ccall "auxi.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
 
 -- | extraction of a submatrix from a complex matrix
 subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double)
 subMatrixC (r0,c0) (rt,ct) x =
-    reshape ct . asComplex . cdat .
+    reshape ct . asComplex . flatten .
     subMatrixR (r0,2*c0) (rt,2*ct) .
-    reshape (2*cols x) . asReal . cdat $ x
+    reshape (2*cols x) . asReal . flatten $ x
 
 -- | Extracts a submatrix from a matrix.
 subMatrix :: Element a
@@ -299,7 +317,7 @@ subMatrix = subMatrixD
 
 diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
     m <- createMatrix RowMajor n n
-    fun // vec v // mat cdat m // check msg [dat m]
+    fun // vec v // matc m // check msg [cdat m]
     return m -- {tdat = dat m}
 
 -- | diagonal matrix from a real vector
@@ -347,11 +365,11 @@ constant = constantD
 
 -- | obtains the complex conjugate of a complex vector
 conj :: Vector (Complex Double) -> Vector (Complex Double)
-conj v = asComplex $ cdat $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1])
+conj v = asComplex $ flatten $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1])
 
 -- | creates a complex vector from vectors with real and imaginary parts
 toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
-toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i]
+toComplex (r,i) = asComplex $ flatten $ fromColumns [r,i]
 
 -- | the inverse of 'toComplex'
 fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double)
@@ -367,7 +385,7 @@ fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
 fromFile filename (r,c) = do
     charname <- newCString filename
     res <- createMatrix RowMajor r c
-    c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
+    c_gslReadMatrix charname // matc res // check "gslReadMatrix" []
     --free charname  -- TO DO: free the auxiliary CString
     return res
 foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 082e09d..dc86484 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -27,11 +27,12 @@ import Foreign.C.Types
 import Data.Monoid
 
 -- | 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
-                  , ptr  :: Ptr t          -- ^ ordinary pointer to the actual starting address (usually the same)
+data Vector t = V { dim  :: Int              -- ^ number of elements
+                  , fptr :: ForeignPtr t     -- ^ foreign pointer to the memory block
                   }
 
+ptr (V _ fptr) = unsafeForeignPtrToPtr fptr
+
 -- | check the error code and touch foreign ptr of vector arguments (if any)
 check :: String -> [Vector a] -> IO Int -> IO ()
 check msg ls f = do
@@ -63,9 +64,7 @@ createVector :: Storable a => Int -> IO (Vector a)
 createVector n = do
     when (n <= 0) $ error ("trying to createVector of dim "++show n)
     fp <- mallocForeignPtrArray n
-    let p = unsafeForeignPtrToPtr fp
-    --putStrLn ("\n---------> V"++show n)
-    return $ V n fp p
+    return $ V n fp
 
 {- | creates a Vector from a list:
 
@@ -80,7 +79,7 @@ fromList l = unsafePerformIO $ do
     f // vec v // check "fromList" []
     return v
 
-safeRead v f = unsafePerformIO $ withForeignPtr (fptr v) $ const $ f (ptr v)
+safeRead v = unsafePerformIO . withForeignPtr (fptr v)
 
 {- | extracts the Vector elements to a list
 
@@ -115,19 +114,14 @@ subVector :: Storable t => Int       -- ^ index of the starting element
                         -> Int       -- ^ number of elements to extract
                         -> Vector t  -- ^ source
                         -> Vector t  -- ^ result
-subVector k l (v@V {dim=n, ptr=p, fptr=fp})
+subVector k l (v@V {dim=n})
     | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
     | otherwise = unsafePerformIO $ do
         r <- createVector l
-        let f = copyArray (ptr r) (advancePtr p k) l >> return 0
-        f // check "subVector" [v]
+        let f = copyArray (ptr r) (advancePtr (ptr v) k) l >> return 0
+        f // check "subVector" [v,r]
         return r
 
-subVector' k l (v@V {dim=n, ptr=p, fptr=fp})
-    | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
-    | otherwise = v {dim=l, ptr=advancePtr p k}
-
-
 {- | Reads a vector position:
 
 @> fromList [0..9] \@\> 7
@@ -149,23 +143,23 @@ join :: Storable t => [Vector t] -> Vector t
 join [] = error "joining zero vectors"
 join as = unsafePerformIO $ do
     let tot = sum (map dim as)
-    r@V {fptr = p, ptr = p'} <- createVector tot
+    r@V {fptr = p} <- createVector tot
     withForeignPtr p $ \_ ->
-        joiner as tot p'
+        joiner as tot (ptr r)
     return r
   where joiner [] _ _ = return ()
-        joiner (V {dim = n, fptr = b, ptr = q} : cs) _ p = do
-            withForeignPtr b  $ \_ -> copyArray p q n
+        joiner (r@V {dim = n, fptr = b} : cs) _ p = do
+            withForeignPtr b  $ \_ -> copyArray p (ptr r) n
             joiner cs 0 (advancePtr p n)
 
 
 -- | transforms a complex vector into a real vector with alternating real and imaginary parts 
 asReal :: Vector (Complex Double) -> Vector Double
-asReal v = V { dim = 2*dim v, fptr =  castForeignPtr (fptr v), ptr = castPtr (ptr v) }
+asReal v = V { dim = 2*dim v, fptr =  castForeignPtr (fptr v) }
 
 -- | transforms a real vector into a complex vector with alternating real and imaginary parts
 asComplex :: Vector Double -> Vector (Complex Double)
-asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v), ptr = castPtr (ptr v) }
+asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }
 
 ----------------------------------------------------------------
 
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index b53d9b7..7f83bcf 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -125,11 +125,26 @@ int multiplyR(int ta, KRMAT(a), int tb, KRMAT(b),RMAT(r)) {
     KDMVIEW(a);
     KDMVIEW(b);
     DMVIEW(r);
+    int k;
+    for(k=0;k<rr*rc;k++) rp[k]=0;
+    int debug = 0;
+    if(debug) {
+        printf("---------------------------\n");
+        printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
+        printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
+        printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
+    }
     int res = gsl_blas_dgemm(
          ta?CblasTrans:CblasNoTrans,
          tb?CblasTrans:CblasNoTrans,
          1.0, M(a), M(b),
          0.0, M(r));
+    if(debug) {
+        printf("--------------\n");
+        printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
+        printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
+        printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
+    }
     CHECK(res,res);
     OK
 }
@@ -140,14 +155,36 @@ int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r)) {
     KCMVIEW(a);
     KCMVIEW(b);
     CMVIEW(r);
+    int k;
     gsl_complex alpha, beta;
     GSL_SET_COMPLEX(&alpha,1.,0.);
     GSL_SET_COMPLEX(&beta,0.,0.);
+    //double *TEMP = (double*)malloc(rr*rc*2*sizeof(double));
+    //gsl_matrix_complex_view T = gsl_matrix_complex_view_array(TEMP,rr,rc);
+    for(k=0;k<rr*rc;k++) rp[k]=beta;
+    //for(k=0;k<2*rr*rc;k++) TEMP[k]=0;
+    int debug = 0;
+    if(debug) {
+        printf("---------------------------\n");
+        printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
+        printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
+        printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
+        //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
+    }
     int res = gsl_blas_zgemm(
          ta?CblasTrans:CblasNoTrans,
          tb?CblasTrans:CblasNoTrans,
          alpha, M(a), M(b),
-         beta, M(r));
+         beta, M(r)); 
+         //&T.matrix);
+    //memcpy(rp,TEMP,2*rr*rc*sizeof(double));
+    if(debug) {
+        printf("--------------\n");
+        printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
+        printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
+        printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
+        //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
+    }
     CHECK(res,res);
     OK
 }