subMatrix

4 files changed, 68 insertions, 11 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index a6f7f0c..a2a70dd 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -98,6 +98,8 @@ createMatrix order r c = do
     p <- createVector (r*c)
     return (matrixFromVector order c p)
 
+reshape c v = matrixFromVector RowMajor c v
+
 transdataG :: Storable a => Int -> Vector a -> Int -> Vector a 
 transdataG c1 d c2 = fromList . concat . transpose . partit c1 . toList $ d
 
@@ -184,8 +186,41 @@ multiply ColumnMajor a b = trans $ multiplyT ColumnMajor a b
 
 multiplyT order a b = multiplyD order (trans b) (trans a)
 
-multiplyD order a b 
+multiplyD order a b
     | isReal (baseOf.dat) a = scast $ multiplyAux order cmultiplyR (scast a) (scast b)
     | isComp (baseOf.dat) a = scast $ multiplyAux order cmultiplyC (scast a) (scast b)
     | otherwise             = multiplyG a b
 
+----------------------------------------------------------------------
+
+-- | extraction of a submatrix of a real matrix
+subMatrixR :: (Int,Int) -- ^ (r0,c0) starting position 
+           -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
+           -> Matrix Double -> Matrix Double
+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 cdat r // check "subMatrixR" [dat r]
+    return r
+foreign import ccall "aux.h submatrixR"
+    c_submatrixR :: Int -> Int -> Int -> Int -> Double ::> Double ::> IO Int
+
+-- | extraction of a submatrix of a complex matrix
+subMatrixC :: (Int,Int) -- ^ (r0,c0) starting position
+           -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
+           -> Matrix (Complex Double) -> Matrix (Complex Double)
+subMatrixC (r0,c0) (rt,ct) x =
+    reshape ct . asComplex . cdat .
+    subMatrixR (r0,2*c0) (rt,2*ct) .
+    reshape (2*cols x) . asReal . cdat $ x
+
+subMatrix :: (Field a) 
+          => (Int,Int) -- ^ (r0,c0) starting position 
+          -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
+          -> Matrix a -> Matrix a
+subMatrix st sz m
+    | isReal (baseOf.dat) m = scast $ subMatrixR st sz (scast m)
+    | isComp (baseOf.dat) m = scast $ subMatrixC st sz (scast m)
+    | otherwise             = subMatrixG st sz m
+
+subMatrixG (r0,c0) (rt,ct) x = reshape ct $ fromList $ concat $ map (subList c0 ct) (subList r0 rt (toLists x))
+    where subList s n = take n . drop s
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 7dcefeb..6ed9339 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -119,23 +119,30 @@ subVector' k l (v@V {dim=n, ptr=p, fptr=fp})
     | otherwise = v {dim=l, ptr=advancePtr p k}
 
 
-{-
+
 -- | creates a new Vector by joining a list of Vectors
 join :: Field t => [Vector t] -> Vector t
-join [] = error "joining an empty list"
+join [] = error "joining zero vectors"
 join as = unsafePerformIO $ do
-    let tot = sum (map size as)
-    p <- mallocForeignPtrArray tot
-    withForeignPtr p $ \p ->
-        joiner as tot p
-    return (V tot p)
+    let tot = sum (map dim as)
+    r@V {fptr = p, ptr = p'} <- createVector tot
+    withForeignPtr p $ \_ ->
+        joiner as tot p'
+    return r
   where joiner [] _ _ = return ()
-        joiner (V n b : cs) _ p = do
-            withForeignPtr b  $ \b' -> copyArray p b' n
+        joiner (V {dim = n, fptr = b, ptr = q} : cs) _ p = do
+            withForeignPtr b  $ \_ -> copyArray p q 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) }
+
+-- | 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) }
+
 constantG n x = fromList (replicate n x)
 
 constantR :: Int -> Double -> Vector Double
diff --git a/lib/Data/Packed/Internal/aux.c b/lib/Data/Packed/Internal/aux.c
index da36035..01a2bb3 100644
--- a/lib/Data/Packed/Internal/aux.c
+++ b/lib/Data/Packed/Internal/aux.c
@@ -89,6 +89,19 @@ int transC(KCMAT(x),CMAT(t)) {
 }
 
 
+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 &&
+            rr==r2-r1+1 && rc==c2-c1+1,BAD_SIZE);
+    DEBUGMSG("submatrixR");
+    KDMVIEW(x);
+    DMVIEW(r);
+    gsl_matrix_const_view S = gsl_matrix_const_submatrix(M(x),r1,c1,rr,rc);
+    int res = gsl_matrix_memcpy(M(r),M(S));
+    CHECK(res,res);
+    OK
+}
+
+
 int constantR(double * pval, RVEC(r)) {
     DEBUGMSG("constantR")
     int k;
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/aux.h
index f45b55a..59d90db 100644
--- a/lib/Data/Packed/Internal/aux.h
+++ b/lib/Data/Packed/Internal/aux.h
@@ -19,3 +19,5 @@ int constantC(gsl_complex *val, CVEC(r));
 
 int multiplyR(int ta, KRMAT(a), int tb,  KRMAT(b),RMAT(r));
 int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r));
+
+int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));